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SECTION  1 


INTRODUCTION 

The  benefits  and  shortcomings  of  micromechanical  analyses  of  uni¬ 
directional  composites  have  now  been  weighed  for  over  a  decade.  A 
major  restriction  has  been  the  problem  of  relating  microanalyses  to 
the  overall  behavior  of  a  laminate.  To  gain  insight  into  this  prob¬ 
lem,  it  is  instructive  to  consider  the  well-known  process  of  design¬ 
ing  with  metals  as  opposed  to  the  task  of  employing  advanced  compo¬ 
sites  in  structures. 

For  centuries  metals  have  been  the  materials  considered  almost 
universally  for  countless  applications.  Metals  with  the  desired  stiff 
ness,  strength,  toughness,  etc.,  can  be  chosen  from  information  at 
arm's  reach,  for  each  specific  application.  This  vast  amount  of  in¬ 
formation  has  greatlv  simplified  design  using  metals.  But  because 
metals  are  in  most  cases  isotropic  and  ultimate  strengths  pertain  to 
all  material  directions,  there  is  often  unneeded  strength  in  certain 
material  directions.  A  truss  member  to  be  loaded  axially  in  tension, 
for  example,  requires  high  strength  in  the  axial  direction,  but  trans¬ 
versely  only  minimal  strength  is  needed.  Thus,  the  high  degree  of  ani 
sotropv  typically  exhibited  by  a  composite  material  is  not  necessarily 
a  handicap  in  many  high-performance  structures.  Composite  materials 
can  be  chosen  for  weight  savings,  and  designed  to  specific  strength 
requirements  in  specific  material  directions.  Because  of  the  vast 
possibilities  created  through  material  design  for  specific 
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applications ,  the  use  of  composite  materials  has  initiated  a  new  era  of 
materials  technology,  and  the  result  is  that  a  much  more  complete  ap¬ 
proach  to  structural  mechanics  has  become  necessary. 

The  transition  from  the  use  of  simple  isotropic  metals  to  ortho¬ 
tropic  composite  plies  and  anisotropic  laminates  for  structural  sys¬ 
tems  necessitates  much  more  involved  analysis  methods.  While  design 
with  metals  has  been  practical  for  centuries,  composite  materials  are 
just  emerging  from  their  own  "Iran  Age."  A  large  amount  of  analysis 
and  technology  must  Le  developed  before  the  full  potential  of  compo¬ 
site  materials  can  be  realized. 

While  overall  properties  of  a  material  are  needed  for  the  design 
of  any  part,  an  additional  dimension  of  composite  materials  behavior 
is  the  unique  problem  of  the  complex  stress  state  present  on  the  micro¬ 
scopic  scale.  When  a  fiber  of  high  strength  and  modulus  is  imbedded 
in  a  matrix  material  of  relatively  low  strength  and  modulus,  local 
stress  concentrations  are  induced  that  are  significant  to  the  overall 
behavior  of  the  composite.  The  mathematical  character ication  of  this 
situation  is  known  as  inicromeehanics  analysis  while  a  larger  scale  an- 
analysis,  i.e.,  lamination  theory,  is  known  as  macromec.hanics . 

Macromechanics  is  concerned  with  the  overall  properties  of  a 
composite  laminate,  or  can  be  concerned  with  the  contribution  of  each 
ply  to  the  overall  effect.  In  contrast,  micromechanics  attempts  to 
predict  such  phenomena  as  matrix  yielding,  crack  initiation,  and  uni¬ 
directional  composite  pl;.r  behavior,  using  an  analvsis  concerned  with 
the  local  scale.  This  is  important  for  understanding  particular 
failure  initiation  modes  and  predicting  properties  of  unidirectional 
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composites.  Most  existing  micromechanics  analyses  employ  a  finite 
element  numerical  scheme,  using  the  material  properties  of  fiber  and 
matrix  to  predict  the  microstress  distributions,  assuming  a  state  of 
generalized  plane  strain.  This  type  of  analysis  could  also  be  performed 
using  a  three-dimensional  finite  element  analysis,  although  much  more 
computer  time  and  memory  would  be  necessary.  Ideally,  a  three-dimen¬ 
sional  micromechanics  analysis  should  eventually  be  developed  as  a 
check  of  the  generalized  plane  strain  formulation.  While  previously 
existing  analyses  allow  only  a  longitudinal  normal  load  in  the  out-of¬ 
plane  direction,  the  present  analysis  permits  a  longitudinal  shear 
loading  capability.  Miller  and  Adams  [1] predicted  in  1977,  "more 
attempts  will  be  forthcoming  to  marry  micromechanical  analyses  .  .  . 
with  various  Lamination  analyses,"  and  this  is  one  of  the  primarv 
reasons  for  the  addition  of  a  longitudinal  shear  loading  capability. 
Results  of  a  Laminate  point  stress  analysis  can  be  analyzed  further 
by  using  such  a  microraechanics  formulation.  The  Laminate  analvsis 
will  reveal  the  stress  state  in  each  ply,  from  which  the  microanalysis 
can  predict  local  inelastic  behavior  and  reveal  the  complex  stress 
state  within  the  matrix  material.  This  type  of  analysis  can  be  con¬ 
sidered  to  be  an  important  first  step  in  bridging  the  gap  between 


raicromechanies  and  macromechanics. 


SECTION  2 


HISTORICAL  REVIEW 

The  importance  of  longitudinal  shear  stress  in  a  composite  material 
is  evidenced  by  the  fact  that  load  is  transferred  to  a  fiber  predomi¬ 
nantly  through  longitudinal  shear  loading.  This  stress  also  happens  to 
act  in  a  weak  direction  of  the  composite,  making  it  a  critical  load. 

Even  though  longitudinal  shear  loading  is  an  important  consideration, 
there  have  been  few  studies  of  it  pertaining  to  micromechanics  of  a 
composite. 

Early  in  1967,  Adams  and  Doner  [2]  revealed  a  numerical  formulation 
involving  the  theory  of  elasticity.  They  modeled  one  quadrant  of  a 
repeating  fundamental  region  of  a  rectangular  array  of  fibers  by  employ¬ 
ing  a  finite  difference  representation,  and  subsequently  solved  the 
on.  blem  by  an  >ver -relaxation  procedure.  Stress  concentration  factors 
and  composite  -shear  moduli  were  calculated  for  various  fiber  volumes 
for  a  number  of  cross-sectional  shapes  of  fibers  in  an  epoxy  matrix. 

This  first  step  was  soon  multiplied  as  numerical  analyses  developed 
and  computers  became  more  advanced. 

A  few  years  later,  a  closed  form  series  solution  was  developed  by 
Sendeckyj  [3]  for  longitudinal  shear  loading.  Admittedly,  the  solu¬ 
tion  was  tedious  due  to  the  solution  technique  employed,  and  fell 
short  of  being  exact  due  to  the  required  truncation  of  the  infinite 
series.  Nonuniform  fiber  spacing,  various  filament  radii  and 
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variation  of  the  shear  modulus  from  fiber  to  fiber  were  some  of  the 
impressive  capabilities. 

Another  elastic  solution  was  achieved  by  Foye  [4]  in  1968.  The 
finite  element  numerical  method  was  employed  in  this  sweeping  genera¬ 
lized  plane  strain  study,  which  included  two  fiber  arrangements,  sepa¬ 
rate  and  combined  loading  of  five  of  the  six  components  of  stress,  con¬ 
tours  of  stresses  in  the  matrix  around  a  fiber,  unidirectional  ply 
composite  properties,  and  an  evaluation  of  the  accuracy  of  various 
finite  element  models.  In  addition,  an  incremental  inelastic  analysis 
was  proposed,  which  was  eventually  employed  by  Baker  and  Foye  [5]  in 
19b9.  This  extended  work  revealed  a  more  legitimate  stress  distribu¬ 
tion  because  inelastic  behavior  was  considered.  Foye  continued  to 
publish  results  of  this  analysis  in  1970  [6]  and  1973  [7],  The  work 
documented  in  1973  essentially  clarified  the  aforementioned  work  of 
1969  [3] . 

Although  this  analysis  was  a  significant  achievement,  there  were 
still  some  limitations  to  be  overcome.  The  iterative  scheme  inherently 
accumulated  error  during  inelastic  increments  which  would  grow  to 
significant  size  as  the  number  of  inelastic  increments  increased.  The 
iterative  inelastic  analysis,  termed  the  method  of  initial  strains, 
had  been  chosen,  even  though  it  degenerates  for  highly  inelastic  beha¬ 
vior,  because  the  alternative  method  (the  tangent  modulus  method) 
would  have  required  an  unavailable  amount  of  computer  memory  for  the 
additional  longitudinal  shear  loading  capability.  Even  though  the 
method  of  initial  strains  had  the  advantage  of  requiring  only  one 
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initial  inversion  of  the  stiffness  matrix,  the  tangent  modulus  method  was 
found  to  be  slightly  faster,  with  equal  accuracy  and  the  ability  to  model 
highly  inelastic  materials  [8]. 

In  the  past  15  years,  computers  have  been  developed  considerably, 
and  the  disadvantages  attributed  to  the  tangent  modulus  method  have  been 
blunted  by  the  increased  size  of  computer  memory.  Advantages  of  the  tan¬ 
gent  modulus  approach  began  to  clearly  emerge  in  the  analysis  presented 
by  Adams  [9,  10]  in  1970.  This  was  further  developed  in  his  subsequent 
work  reported  in  References  [11,  12,  13,  14).  This  analysis  method  was 
subsequently  adapted  by  Miller  and  Adams  [1,  15,  16],  incorporating  work 
by  Crossman  [17]  and  Sranca  [18].  It  was  a  generalized  plane  strain 
formulation,  including  longitudinal  normal  loading,  but  not  longitudinal 
shear.  The  work  of  Miller  and  Adams  was  particularly  valuable  due  to 
the  addition  of  a  hygrothermal  loading  capability  along  with  hygrother- 
mally-dependent  material  properties,  and  3ranca's  [18]  efficient  loading 
scheme.  This  generalized  plane  strain  approach  is  readily  expanded  to 
include  a  longitudinal  shear  loading  capability. 


SECTION  3 


ANALYSIS  METHOD 

The  changes  in  theory  that  arise  due  to  the  incorporation  of  an 
out-of-plane  shear  capability  in  a  generalized  plane  strain  formula¬ 
tion,  which  includes  only  a  normal  stress  in  the  out-of-plane  direc¬ 
tion,  are  significant.  In  addition,  there  are  far  reaching  conse¬ 
quences  in  the  finite  element  solution  technique.  The  added  stress 
and  strain  components  require  revisions  in  the  elastoplastic  consti¬ 
tutive  equations.  In  addition,  the  special  treatment  of  the  compati¬ 
bility  of  additional  boundary  conditions  results  in  a  special  strain- 
displacement  reLation.  Implementing  these  changes  into  a  finite 
element  analvsis,  although  following  classical  developments,  is 
difficult.  Use  of  the  existing  computer  program  [1]  required  its 
thorough  revision.  In  many  cases,  where  before  a  plane  stress  state 
was  the  only  consideration,  a  stress  tensor  was  involved  in  developing 
suitable  failure  criteria,  principal  stress  calculations,  etc. 

A  complete  description  of  each  aforementioned  consideration  is 
presented  in  the  following  subsections. 

Generalized  Plane  Strain 

Leknitskii  [ Id]  defines  generalized  plane  strain  in  a  very  general 
manner,  vhich  has  been  simplified  for  purposes  of  the  present  analysis. 
This  treatment  allows  displacements  to  occur  in  all  three  coordinate 
iireetions,  yet  retains  the  advantages  of  the  plane  strain  assumption. 
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Specifically,  each  displacement  is  dependent  upon  the  x-  and  '/-coor¬ 
dinate  directions,  and  the  displacement  in  the  z-direction  has  an 
added  linear  dependence  on  the  z-coordinate ,  which  is  considered  the 
axial  coordinate  of  a  composite  in  the  present  analysis  (Figure  1). 


Figure  1.  Fiber  packing  arrangement  of  a  unidirectional  composite. 

The  displacement  functionals  in  equation  form  are: 

u  =  u  ( x ,  y ) 

v  =  v(x,y)  (1) 

w  =  w  v  x ,  y )  +  C ;  z 

where 

u  represents  the  x-displacement 
v  represents  the  y-displacement 
w  represents  the  z-displacement 
and  Ci  is  a  yet  unspecified  constant. 
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Previous  investigators  [1]  making  use  of  the  generalized  plane  strain 
assumption  simplified  the  expression  for  the  z-displacement  by  assuming 
it  was  dependent  only  on  the  axial  coordinate  position.  Eliminating  the 
functional  dependence  on  the  x-  and  y-coordinates  essentially  eliminated 
axial  shear  deformations,  while  allowing  only  constant  axial  normal  dis¬ 
placements. 

Including  the  x  and  v  dependence  of  the  z-displacement  allows  a 
special  form  of  axial  shear  deformation  corresponding  to  the  generalized 
plane  strain  treatment.  Expressions  for  strain  can  now  be  calculated 
for  the  continuum,  and  simplified  according  to  Eqs.  (1): 


3u  _  £U  _3v 

x  3x  Yxy  3v  3x 


where  t  represents  normal  strain  and  y  represents  engineering  shear 
strain.  These  expressions  govern  the  displacement  of  the  continuum  in 
question. 

Constitutive  Equations 

The  key  element  of  the  relationship  between  load  and  displacement 
of  a  continuum  is  the  consititutive  equation.  Each  constituent  material 
of  a  composite  has  unique  properties  represented  in  its  own  constitu¬ 
tive  equation.  In  the  present  analysis,  the  fiber  material  is  consid¬ 
ered  to  be  isotropic  and  transversely  isotropic  in  order  to  model 
fibers  such  as  graphite,  but  can  be  reduced  to  a  simple  elastic, 
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isotropic  material  for  the  purpose  of  characterizing  fibers  such  as 
boron  and  glass.  The  matrix  material  is  considered  to  be  isotropic 
and  elastoplastic ,  the  plastic  response  being  modeled  by  the  Prandtl- 
Reuss  flow  rule. 

A  constitutive  equation  involving  only  four  stress  and  strain 
components  has  been  derived  in  detail  [1],  but  the  two  additional 
longitudinal  shear  components  require  additional  consideration .  As 
discussed  by  Baker  and  Foye  [5],  the  two  additional  shear  stress- 
shear  strain  equations  are 


’xz  =  r/ 

yxz 

'  v  z  —  G f 

7  V  2 

(3) 

where  G'  is  the  longitudinal  shear  modulus  and  t  represents  shear 
stress.  When  Eqs.  (3)  are  included  in  the  complete  set  of  constitutive 
equations  for  transversely  isotropic  elastic  behavior,  a  material 
properties  matrix  (Di  is  generated  as 
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•flow  rule.  Reference  [ !  |  discusses  sources  of  more  information  related 
to  these  two  criteria.  A  plane  strain  treatment  of  elastoplastic 
b,  savior  by  Adams  [11  j  is  readily  adapted  to  a  generalized  plane  strain 
treatment ,  as  outlined  irieiiv  in  Appendix  A.  The  resulting  constitu- 

r.  Lve  caation  is 


.c  :  ud  :  v  idti.il  :  ..  . .  <  scribed  m  Appendix  A. 


A  •<a;..r  oousi.h-rut  ion  in.  cue  finite  element  formulation  is  the 
•  I.-  v  e  a  combine:  .  nil  at  capability  ,  whore  anv  combination 

load  can  bo  .vpn I  ;«.d  jprr.u !  t atioouslv .  This  leads  to  the  seemingly 
i;r>ns-<  ip. le  to-.'.:  of  1  ins  ..-•narare  and  contradictory  boundary  con- 
.  :  ton--.  *i mult ancon:  !••.  Before  these  factors  can  be  considered  in 
i :  !  ,  •  :«•  r'ur.d  amen  i.  a  1  problem  must  first  be  set  up. 

IS.e  finite  element  model  to  no  used  in  this  investigation  is 
pi  t;:e  ussunpt  ion  t ha t  the  fibers  are  distributed  in  a  square 
irr.iv ,  !.  d-iowti  in  Figure  1,  although  anv  rectangular  array  is  entirely 

: i . i . ■ .  previous  invest  •;  gators  [1,  !<i .  20)  have  shown  that  there  is 
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little  loss  of  modeling  accuracy  due  to  this  assumption.  Considering 
the  behavior  of  a  unidirectional  composite  with  a  relatively  long  axial 
dimension,  certain  load  symmetry  and  geometric  symmetry  assumptions  can 
be  made,  as  described  in  Reference  [4].  While  there  are  four  axes  of 
geometric  symmetrv  for  the  unit  cell  (Figure  2a),  only  two  can  be  used 
due  to  load  symmetry  considerations .  When  considering  that  each  load 
is  assumed  to  deform  the  unit  cell  in  a  uniform  manner,  load  symmetry 
about  the  x  and  y  axes  for  all  five  stress  components  can  be  easily 
seen.  Only  one  quadrant  of  the  unit  cell  (Figure  2b)  need  be  con¬ 
sidered  to  describe  the  behavior  of  the  unit  cell  and  of  an  entire 
continuum  of  unit  cells. 

A  detailed  description  of  the  boundary  conditions  is  necessary, 
beginning  with  the  normal  displacements.  Due  to  the  previously 
mentioned  constraints  of  the  system,  normal  displacements  of  the 
boundaries  of  tne  quadrant  (Figure  2b)  are  restricted  to  those  which 
cause  the  boundarv  to  displace  only  parallel  to  the  original  boundary, 
ihis  being  explained  in  great  detail  elsewhere  [1],  attention  will  be 
focused  on  the  shear  displacements.  In  the  absence  of  longitudinal 
shear  loading,  deformation  in  the  z-direction  is  simply  a  constant, 
uniform  deformation  of  the  entire  frontal  cross-section  cf  the  quadrant 
Shear  deformations  constrain  only  the  loaded  boundary  to  displace  uni¬ 
formly  in  the  z-direction,  while  the  opposite  face  remains  stationary 
in  the  z-direction.  A  r  shear  stress,  for  example,  causes  the  face 
at  v  =  b  'see  Figure  S')  to  move  uniformly  in  the  z-direction  while  the 
face  at  y  =  0  is  fixed  in  the  z-direction.  Meanwhile,  the  faces  at 


a) 


\r 


b)  Quadrant  to  be  analvzed 


Figure 


Symr.etr”  simplification  of  the  fiber  arrangement 
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Figure  3.  Shear  deformation  of  the  finite  element  model  due  to 
shear  loading. 


x  =  0  and  x  =  a  are  free  to  move  in  the  z-direction  as  strains  are 
induced.  Similar  treatment  of  the  other  longitudinal  shear  stress 
follows.  It  will  be  noted  that  these  restrictions  are  not  compatible, 
i.e.,  while  the  face  at  x  =  0  is  required  to  be  fixed  in  one  case,  it 
is  required  to  be  free  in  the  other  case.  To  solve  this  anomaly,  the 
analysis  must  be  considered  for  the  representative  elements  of  the 
entire  quadrant. 

Constraints  following  those  of  Reference  [2]  are  imposed  on  the 
special  generalized  plane  strain  deformation  of  the  finite  element 

y  *  b 
y  *  0 

x  »  0,  x  ■  a 


continuum  due  to  a  tvz  shear  stress: 

w  =»  C  along 

w  =  0  along 
^  w 

^  =  0  along 


(7) 
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Likewise  for  txz  loading: 


w  =  C 

along  x  = 

a 

w  =  0 

along  x  = 

0 

(8) 

|S-  0 

3y 

along  y  * 

0,  v  =  b 

where  C  represents  a  constant  displacement .  It  will  again  be  noted 
that  the  two  shear  boundary  conditions  cannot  be  applied  to  the  same 
node  simultaneously.  A  method  of  considering  the  two  shear  loads  in 
separate  problems  while  permitting  combined  shear  loading  is  desired 
in  the  finite  element  formulation. 

It  is  possible  to  apply  the  two  shear  boundary  conditions  separ¬ 
ately  by  adjusting  the  elemental  strain-displacement  relation.  The 
derivation  of  this  relation  for  a  regular  treatment  of  generalized 
plane  strain  is  outlined  in  Appendix  B.  The  resulting  Eq .  (B-17)  must 
be  adapted  to  manage  the  shear  boundary  conditions  separately.  Since 
the  boundary  conditions  are  concerned  with  displacements,  the  shear 
displacetiR  nts  must  be  considered  to  be  separate  and  unique  entities  at 
each  node  point  rather  than  combined  into  a  general  z-di splacement  as 
in  Eq .  (B-17).  Thus,  a  t  shear  stress  induces  a  z-disp lacement  w  , 
and  a  shear  stress  induces  a  z-displacement  wv,-  These  separate 
displacements  can  be  arbitrarily  specified  without  affecting  other 
boundary  conditions  if  the  strain-displacement  matrix  equation  is 
expanded  slightly  as  follows: 
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(9) 


where  i,  j,  and  k  represent  the  three  respective  nodes  of  an  ele¬ 


ment,  and  w  represents  the  constant  normal  axial  displacement. 

It  will  be  noted  that  the  three  extra  columns  in  the  matrix  and  the 
added  zeros  effectively  prevent  the  longitudinal  shear  strains  from 
influencing  each  other  and  the  remainder  of  the  strains.  The  added 
terms  in  the  displacement  vector  will  also  be  noted.  Observing  the 
expression  for  elemental  shear  strain  yy2  above,  provides  a  check  of 
the  last  boundary  condition  in  Eq.  (7)  which  is  now  inherent  in  the 

Cf  W 

strain-displacement  equation.  The  shear  strain  YyZ,  equal  to  —  ,  is 

seen  to  depend  only  on  the  variables  c^,  and  c^,  which  are  merely 

differences  in  the  x-coordinates  of  the  nodes  of  an  element.  Because 

the  w  displacement  is  treated  separately  from  the  w  displacement, 
xz  yz 


the  expression  for  y  does  not  include  a  y 

v  VZ  xz 


dw 


rxi 


shear  strain,  and 


this  is  inherent  in  the  expression.  Therefore,  this  special  boundary 
condition  has  been  incorporated  into  the  strain-displacement  equation. 


t 
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The  remaining  finite  element  formulation  follows  a  tnree- 
dimensional  formulation  by  Zienkiewicz  [21],  as  shown  in  abbreviated 
form  in  the  following  equations: 


[k]e  =  JvoltB]J  [D][B]  d(vol) 

(10) 

{ f  }  =  (Kp{F} 

(ID 

U)  =  [B]  {6} 

(12) 

i  J s  =  [D][B]  { ij 

(13) 

where 

{6;  =  nodal  displacements 

(c)  =  elemental  strains 

[ B ]  =  strain-displacement  matrix 

[ D ]  =  constitutive  equation 

•.a;  =  elemental  stresses 

[K]e  =  elemental  stiffness  matrix 

{F;  =  nodal  force  vector 

[K]  =  inverse  of  the  global  sciffness  matrix. 

Elemental  stiffness  matrices  are  developed  in  Eq.  i,10)  ,  which  combined 
together  for  the  entire  model  form  a  global  stiffness  matrix.  This 
is  inverted  Cor  use  in  Eq .  (11),  to  solve  for  nodal  displacements. 
Element  strains  and  stresses  are  calculated  from  the  displacement 
vector  in  Eqs.  (12)  and  (13),  respectively.  Two  differences  from 
Zienkiewicz ' s  formulation  are  the  nodal  force  vector  and,  of  course, 
the  nodal  displacement  vector.  The  displacement  vector  was  previously 
defined;  the  nodal  force  vector  proves  to  be  very  similar,  due  to  the 
fact  that  there  are  five  force  components  and  displacement  components 
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possible  at  each  node  of  the  finite  element  grid.  Because  the  solution 
technique  involves  the  force  vector  to  a  large  extent,  it  will  be 
described  in  more  detail  next. 

Solution  Technique 

A  finite  element  analysis  computer  program  developed  by  Miller 
and  Adams  [1]  has  been  adapted  for  the  modified  form  of  generalized 
plane  strain.  The  additional  components  of  stress  and  strain  result 
in  the  necessity  of  a  thorough  revision  of  the  basic  solution  technique. 
The  Branca  technique  [IS]  for  applying  loads  and  boundary  conditions, 
and  the  specialized  Gaussian  elimination  solution  of  the  stiffness 
matrix,  are  the  major  burdens  in  the  revision.  Developing  these  new 
techniques  revealed  a  more  efficient  procedure  for  storing  the  stiff¬ 
ness  matrix.  The  solution  technique  involves  the  application  of  loads 
and  boundary  conditions,  which  is  a  logical  starting  point  for  the 
following  description. 

To  describe  the  Branca  technique  [18],  a  small  two-element  model 
shown  in  Figure  4  will  be  considered.  The  applied  loads  shown  are 


Figure  4.  Simplified  model  for  example  solution  technique. 
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actually  applied  as  nodal  forces,  resulting  in  five  separate  loading 
components  possible  at  each  node.  As  explained  in  detail  by  Miller 
and  Adams  [1],  the  axial  normal  load  is  assumed  to  cause  a  constant 
displacement  in  the  z-direction  at  every  node.  The  solution  technique 
for  that  z-displacement  combines  all  the  axial  stiffness  terms  and 
axial  load  terms  into  one,  and  one  displacement  term  is  solved  for, 
which  represents  the  z-displacement  of  ail  node  points.  Therefore, 
with  the  remaining  four  loads  at  each  node,  the  resulting  load  vector 
for  the  entire  model  is  17  terms  long,  as  is  the  displacement  vector 
(see  Eq.  14),  (where  F i 7  and  *17  represent  the  z-direction  'orce  and 
disolaeement,  respectively). 

The  global  stiffness  matrix  is  created  from  the  elemental  stiff¬ 
ness  matrices  by  means  of  classical  techniques  [21,  22],  but  the 
application  of  boundary  conditions  to  the  global  stiffness  matrix  is 
necessary.  Assume  the  model  in  Figure  4  to  be  loaded  and  constrained 
in  a  manner  similar  to  the  total  finite  element  grid.  Before  boundary 
conditions  are  applied,  the  diagonalized  stiffness  matrix  appears  as 
in  Eq.  (14).  Each  row  of  the  original  matrix  governs  the  behavior  of 
a  node  due  to  a  particular  load.  There  are  four  rows  for  every  node 
in  the  finite  element  grid.  In  addition,  the  last  row  and  column 
represent  behavior  due  to  normal  axial  loads.  It  will  be  noted  that 
terms  in  a  square  stiffness  matrix  are  arranged  in  regular  rows  and 
columns,  but  upon  diagonalization  the  columns  are  skewed  upward  to 
the  right.  A  row  in  a  square  stiffness  matrix  is  represented  differ¬ 
ently  in  a  diagonalized  stiffness  matrix;  terms  that  were  to  the  left 


of  the  diagonal  are  now  in  the  skewed  column,  while  the  remainder  of 
the  row  now  begins  in  Column  1.  A  typical  diagonalized  row  is 
underlined  in  Eq.  (14). 

Table  1  shows  tiie  nodal  constraints  upon  the  model  in  Figure  4 
necessary  to  mimic  the  behavior  of  a  composite  in  a  micromechanics 
analysis.  Effects  upon  the  stiffness  matrix  due  to  a  fixed  node  are 
simple;  the  row  of  the  stiffness  matrix  pertaining  to  the  direction 
of  fixity  of  the  node  is  zeroed  except  for  the  diagonal  term,  which 
is  given  a  value  of  one.  It  will  be  noted  that  part  of  i  zeroed  row 
is  skewed  as  described  previously  and  is  shown  bv  the  underlined  terms 
in  Eq .  (151  .  Upon  solution,  this  causes  a  zero  displacement  to  be 
calculated  at  this  node.  When  all  the  fixed  constraints  are  invoked, 
the  stiffness  matrix  appears  as  in  Eq .  (15). 
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five  boundary  condition  columns  enlarged,  that  were  initially  shown 
in  Eq.  (15).  The  new  terms  were  created  by  columnizing,  for  example; 
Rows  9  and  13  of  Eq .  (15),  representing  the  x-direction  stiffnesses 
of  Node  3,  are  now  stored  in  the  X  Column  of  Eq .  (16).  The  previous 
positions  of  the  terms  that  were  just  relocated  must  be  eliminated 
since  their  effects  have  been  taken  into  account  in  the  boundary 
condition  columns:  setting  diagonal  terms  equal  to  one  and  the 
remaining  terms  in  question  equal  to  zero,  eliminates  their  contri¬ 
bution.  Because  boundary  conditions  were  involved  in  each  direction 
on  each  node,  in  this  case,  the  regular  stiffness  matrix  is  totally 
eliminated  when  its  terms  are  combined  into  the  boundary  condition 
columns,  as  seen  to  the  right  of  the  vertical  dashed  line  (Eq.  16). 

The  second  step  is  shown  by  arrows  in  Eq .  (16),  representing  the 
combination  of  boundary  condition  rows,  located  within  the  special 
boundary  condition  e  mnn  area,  into  the  final  effective  stiffness 
area  (Eq.  17).  1'hese  manipulations  to  include  the  two  additional 
shear  boundary  conditions  being  considered  follow  closely  the  tech¬ 
nique  used  in  References  [18]  and  [1].  For  example,  each  term  of 
Row  6  in  Eq .  (16),  (representing  y-direction  stiffness)  is  added  to 
its  respective  column  in  Row  li  (as  shown  by  the  respective  arrow). 

Eq .  (17)  shows  the  resulting  terms  of  the  summation  due  to  the  combin¬ 
ing  operation.  Final  manipulations  of  the  effective  stiffness  area 
are  employed  to  simplify  the  Gaussian  elimination.  Since  Eq .  (17)  is 
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The  size  of  the  stiffness  matrix  evident  in  Eq.  (14)  becomes 
enormous  when  dealing  with  large  finite  element  models.  This  stiff¬ 
ness  matrix  uses  by  far  the  majority  of  computer  memory  necessary. 

Kid  decreasing  its  size  also  reduces -the  computer  time  needed  to  solve 
it.  One  way  to  decrease  the  size  of  the  stiffness  matrix  is  evident 
in  Eq .  (18),  where  the  large  stiffness  matrix  was  reduced  to  a  system 
of  five  equations.  This  method  eliminates  the  rows  and  columns  of  the 
stiffness  matrix  where  boundary  conditions  are  involved.  These  elimina¬ 
tions  take  place  as  the  stiffness  matrix  is  being  assembled.  The  amount 
of  the  redact i  «  possible  depends  entirely  on  the  number  of  boundary 
conditions  invoked,  out  the  advantage  is  significant  for  the  types  of 
finite  element  cr:  i.,  employee  in  this  investigation,  and  doesn't 
reduce  accuracy .  It  should  be  noted  that  this  method  of  eliminating 
rows  and  columns  was  also  utilized  by  Branca  [18]. 

Next,  general  applied  loads  are  input  to  Eq .  (18)  and  overall 
displacements  arc  •  1  ved  for.  The  resulting  displacements  are  those 
or  tiie  bounder-,  iso,  which  are  then  substituted  back  into  the  original 
displacement  vert  r  i shown  in  Eq .  i4  > .  The  elemental  strains  and 
stresses  can  t.;ien  be  calculated. 

The  sped,;  it  t:  f  ness  formula  tior  discussed  above  caused  com¬ 
plexities  in  the  Gaussian  elimination  procedure.  Branca  [18]  who 
first  developed  the  special  Gaussian  elimination  scheme,  also  testified 
that  it  requir'd  intricate  bookkeeping.  Bookkeeping  complexities  are 
compounded  : taie  p-csent  analysis  due  to  tile  added  boundary  conditions, 
but  the  theory  up.  n  wnic.li  Gaussian  elimination  is  based  remains 
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unchanged.  Thus  the  Gaussian  elimination  procedure  degenerates  into 
an  involved  exercise  in  bookkeeping  and  programming. 

Computer  Implementation 

Implementing  the  preceding  theory  into  the  FORTRAN  computer 
program  of  Reference  [1]  results  in  the  new  generalized  plane  strain, 
finite  element  micromechanics  analysis  described  in  the  flowchart  in 
Appendix  C.  The  incremental  procedure  utilized  by  Miller  and  Adams 
[1]  remains  intact  in  the  present  analysis  except  for  the  calculation 
of  octahedral  shear  stress,  which  now  includes  the  two  added  shear 
stress  terms.  This  tangent  modulus  method  enables  highly  inelastic 
materials  to  be  managed  easily. 

In  developing  the  analysis,  another  consideration  was  found  to 
be  essential  for  an  efficient  computer  program.  Designing  the  finite 
element  mesh  efficiently  has  a  profound  effect  on  the  size  of  the 
stiffness  matrix.  The  highest  difference  in  node  numbers  in  any  one 
element  determines  the  bandwidth  of  the  stiffness  matrix,  i.e., 

Bw  =  (R  +  1)4  +  5  (19) 

where 

BW  =  bandwidth  of  stiffness  matrix 
R  =  highest  node  number  difference 
The  addition  of  5  is  due  to  the  special  boundary  condition  columns. 

The  finite  element  mesh,  besides  the  requirement  to  have  a  fine  mesh 
in  areas  of  high  stress  gradients,  is  also  required  to  have  a  minimum 
R  value.  Additionally,  due  to  the  loading  technique,  the  highest  node 
number  is  also  required  to  be  located  at  the  upper  right-hand  corner 
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or  Che  mesh.  A  model  was  developed  as  shown  in  Figure  5  Chat  fulfilled 
these  requirements,  by  arranging  node  points  in  rows  and  columns 
(although  obscure),  and  numbering  the  nodes  one  row  at  a  time  in 
successive  columns.  There  ’"e  13  rows  and  13  columns  resulting  in 
lb9  nodes  and  288  elements.  The  resulting  R  value  is  14,  making  the 
bandwidth  b5,  as  calculated  from  Eq .  (19). 

In  a  typical  program  run,  utilizing  the  finite  element  model  of 
Figure  5,  the  stiffness  matrix  alone  requires  37,180  words  (60  bits 
per  word)  of  central  memory  on  the  CDC  Cyber  760  computer  system.  This 
amount  is  over  40  percent  of  the  entire  central  memory  necessary,  which 
is  approximately  90,000  words.  An  average  time  used  by  the  central 
processor  to  run  the  program  for  each  increment  is  5.57  seconds.  This 
value  varies  between  5  seconds  and  8  seconds  per  increment,  depending 
on  the  number  of  increments,  for  the  applications  in  this  investigation 
and  most  investigations. 

Increment  sizes  should  be  under  500  psi  for  in-plane  normal 
stresses  and  longitudinal  shear  stresses,  while  axial  normal  stresses 
should  be  about  1000  psi  or  under  for  typical  fibers  of  high  stiffness 
(glass  or  graphite).  During  highly  inelastic  behavior,  increments 
should  be  decreased  accordingly  to  values  on  the  order  of  100  psi. 

Failure  Theories 

There  are  three  failure  theories  necessary  for  a  comprehensive 
prediction  of  failure  in  the  present  micromechanics  analysis.  An 
octahedral  shear  stress  criterion  is  used  to  define  failure  of  the 
matrix.  This  is  a  measure  of  the  distortional  energy  stored  in  the 


Figure  5.  Typical  finite  element  model  of  quadrant  to  be  analyzed. 
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matrix.  The  limitation  of  this  criterion  is  that  it  is  not  suitable 
for  a  hydrostatic  load,  since  there  is  no  distortional  component 
present  in  hydrostatic  loading.  Therefore,  a  hydrostatic  criterion 
is  also  invoked  by  testing  for  a  hydrostatic  stress  state  higher  than 
the  ultimate  value. 

In  prior  analyses  [1]  the  fiber  was  tested  for  failure  by  a 
simple  maximum  stress  criterion.  Since  only  three  normal  stresses 
existed,  a  more  complicated  criterion  was  not  needed.  Because  shear 
stresses  also  exist  in  the  present  analysis,  a  new  criterion  must  be 
invoked.  The  fact  that  graphite  is  an  orthotropic  (transversely  iso¬ 
tropic)  material  increases  the  complexity  of  the  problem.  A  criterion 
specifically  designed  for  orthotropic  materials,  and  which  considers 
t tie  entire  stress  tensor,  is  the  Tsai-Wu  [23]  failure  criterion. 
Experimental  results  show  that  the  Tsai-Wu  criterion  predicts  failure 
far  better  than  a  maximum  stress  or  maximum  strain  criterion  [24]. 

The  Tsai-Wu  tensor  form  is  also  of  a  more  general  character  than  the 
Tsai-Hill  criterion  [2a],  The  only  awkward  characteristic  of  the 

Tsai-Wu  criterion  is  the  P,  term.  This  has  been  discussed  widely  in 

1 

the  literature,  but  Narayanswami  and  Adelman  [25]  show  that  neglecting 
this  term  rarely  causes  an  error  greater  than  10  percent.  This  is  the 
so-called  "Modified  Tsai-Wu"  failure  criterion.  Assuming  tensile  and 
compressive  allowable  stresses  to  be  equal,  the  form  of  the  equation 
becomes : 


Pl°l  +  +  P3o3  +  P4o4  +  p5os  +  P6a6 

+  F1  1°  1  :  +  VZ^  22+  P3  3^  3  2  +  P44°<,2+  P55°5 


+ 


P  2 

66^6 


1 


(20) 


32 


where c  represents  the  actual  stress  components,  and  the  P  terms 
represent  stress  allowables  as  defined  by  Tsai  and  Wu  [23], 
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SECTION  4 


CONSTITUENT  MATERIAL  PROPERTIES 

The  matrix  system  used  for  experimental  verification  pur^.-is  in 
the  present  investigation  was  Hercules  3501-6  epoxy  resin  [26].  u.cause 
longitudinal  shear  loading  is  a  major  consideration  in  the  present 
study,  it  is  appropriate  to  employ  matrix  constituent  material  proper¬ 
ties  derived  from  longitudinal  shear  experimental  data.  Solid  rod 
torsion  test  shear  data  were  available  for  this  matrix  material  [27]. 
This  test  has  been  shown  to  be  a  viable  means  of  determining  shear 
moduli  and  shear  strengths  [28].  Each  test  specimen  was  approximately 
4  inches  long,  with  a  diameter  of  \  inch,  fabricated  in  a  mold  similar 
to  those  shown  in  Reference  [28].  These  tests  were  performed  on  dry 
specimens  and  specimens  sa-irated  with  moisture,  at  three  temperature 
conditions,  viz.,  21°C  (room  temperature)  100°C,  and  160°C. 

Data  from  either  shear  tests  or  tensile,  tests  can  be  readily  con¬ 
verted  to  octahedral  shear  stress-octahedral  shear  strain  expressions 
for  input  to  the  analysis  [29].  Previous  investigators  [ 30 ]  had  already 
calculated  the  octahedral  shear  stress-octahedral  shear  strain  behavior 
of  the  epoxv  matrix  using  uniaxial  tensile  test  data,  as  indicated  in 
Figure  6.  SimiLar  results  reduced  from  solid  rod  torsion  test  data  by 
the  present  investigator  are  shown  in  Figure  7.  The  tensile  test  speci¬ 
mens  were  of  a  "dog-bone"  configuration;  approximate  dimensions  were 
1/10  inch  thick,  4  inch  width,  24  inch  gage  length,  and  5  inch  overall 
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length.  To  use  either  of  these  sets  of  data  in  the  micromechanics 
computer  program,  they  must  be  expressed  in  equation  form,  with 
temperature  and  moisture  as  independent  variables.  First,  a  curve¬ 
fitting  method,  developed  by  Richard  and  Blacklock  [31],  was  used  to 
fit  each  data  curve  (at  each  test  condition).  This  method  defines 
each  curve  in  terms  of  an  initial  slope,  an  asymptotic  stress  value, 
and  a  radius  of  curvature  that  connects  these  two  values.  The  general 
form  of  the  Richard-Blacklock  equation  for  stress  is: 


a 


L 


E  c 

!  E  c  in 


(21) 


where 

E  =  initial  slope  of  stress-strain 
=  the  asymptotic  stress  value 
n  =  inverse  of  the  radius  of  curvature 
£  =  strain 

Each  constant  (E,  c  and  n)  from  each  test  condition,  and  values  for 
ultimate  stress,  are  then  used  in  another  curve-fit  relation  which  ex¬ 
presses  each  value  by  an  equation  dependent  upon  temperature  and 
moisture.  It  is  a  second-order,  least-squares  development  resulting  in 
an  equation  with  six  constants.  For  example,  the  initial  slope  would 
be  expressed  as: 

E  =  C,T-  +  C  tM-  +  C3MT  +  CUM  +  CJ  +  C6  (22) 
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where 

T  *  Temperature  ('  C) 

M  »  Percent  moisture  bv  weight 

Resulting  expressions  for  F.,  n  and  ultimate  stress  are  defined  by 

the  constants  listed  in  Fable  2. 


TABLE  2 

EXPERIMENTALLY-DETERMINED  CONSTANTS  FOR  TEMPERATURE-  AND 
MOISTURE-DEI’ENDENT  OCTAHEDRAL  SHEAR  STRESS-OCTAHEDRAL 
SHEAR  STRAIN  CURVES  AS  OBTAINED  FROM  SOLID  ROD  TORSION 
TEST  DATA  FOR  HERCULES  3501-O  EPOXY  RESIN 


Values 
Dependent 
on  T  and  M 

CONSTANTS 

□ 

c  * 
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C  i  , 
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C6 

1 
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M 

O 

4 

O 

10 

o 
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:o*-!  o 
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n 
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vO 

O 

X 
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0  -3 . 399  x 
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•57  x 

:o'*2.348 

1 

cr 

r+- 

P 

6.315  X 

io';,:.773  x  ;o-i 

'3.-78  x  10* :‘-l . 7*6 x 

io: 

301  x 

10;  1 1 . 322  x 

£1 

Incorporating  these  equations  directly  into  the  computer  program  en¬ 
ables  the  calculation  of  tangent  modulus  vali es  for  each  load  incre¬ 
ment,  as  needed  for  the  inelastic  material  properties  matrix  defined 
in  Section  3.  The  ultimate  octahedral  shear  stress  is  also  calculated 
for  use  in  the  octahedral  shear  stress  failure  criterion,  also  dis¬ 
cussed  in  Section  3, 
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Comparing  Clio  octahedral  shear  stress-octahedral  shear  strain 
curves  generated  from  tensile  and  from  shear  test  data.  Figures  6  and 
7,  respec t ive 1 v,  the  shear  test  data  show  modulus  values  slightly 
higher,  hut  close  enough  to  attest  to  the  general  accuracy  of  the  ex¬ 
perimental  data.  Another  observation  is  the  significant  extent  of  in¬ 
elastic  behavior  of  the  epoxy  resin  when  subjected  to  shear  test  con¬ 
ditions,  (Figure  71  in  contrast  to  the  limited  plastic  response  ex¬ 
hibited  by  the  uniaxial  tensile  data.  A  possible  reason  for  this 
difference  is  the  unique  effect  of  flaws  in  each  type  of  specimen.  The 
thin,  "dog-bone",  tensile  specimen  is  more  susceptible  to  stress  con¬ 
centrations  duo  to  flaws  such  as  voids  or  microcracks  than  the  shear 
specimen.  That  is,  the  shear  specimen  is  more  stable  in  the  presence 
of  flaws  because  only  the  material  on  the  surface  of  the  rod  is  highlv 
stressed.  This  surface  is  smooth  and  has  i  low  number  of  flaws  be¬ 
cause  the  specimen  is  fabricated  in  a  careful iv  polished  steel  mold. 
Internal  f  laws  in  the  rod  sr><  c  i.mens  have  a  reduced  influence  compared 
to  those  m  the  tensile  specimens,  which  are  subjected  to  a  uniformlv 
high  stress  throughout  their  volume. 

Another  f.H'tor  possihl  '  contributing  to  the  difference  is  the 
ic  curacy  at  the  test  eondit  ions ,  r  he  scour. u-v  and  untf’rmity  of 

test  t ..  taper  at  ure  mu  moisture  conditions  which  were  maintained.  In 
both  cases,  mo ist ur *-cond 1 1 1 lined  specimens  were  ac tun ' 1 v  tested  in  a 
dry  environment,  winch  tended  to  dr.ve  cut  the  surface  moisture,  es- 
ci'i  tallv  at  the  high  test  temperature;-.  '.'he  stress  induced,  in  the 
t-.-asi  ie  spe-  mens  due  to  this  mo:  ur.  grn.t ut  ;r  ban  1  v  if  :  octet!  the 
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modulus  and  strength  much  more  than  in  the  thicker  shear  specimens. 

Also,  the  use  of  "clam-shell"  heaters  for  elevated  temperature  tests 
was  a  possible  source  of  error  in  both  test  cases  because  of  the  po¬ 
tential  nonuniformity  of  temperature  around  the  test  specimen,  and  the 
difficulty  in  accurately  controlling  the  temperature  in  such  a  setup. 

It  is  difficult  to  attribute  the  differences  in  behavior  between 
the  tensile  and  shear  specimens  to  experimental  error  alone,  however, 
and  at  this  point  it  is  impossible  to  use  one  form  of  the  data  to  refute 
the  other.  Unfortunately,  there  are  no  other  published  shear  data 
available  that  could  provide  insight  to  the  problem.  However,  there  are 
some  tensile  data  available  that^requi re  a  review. 

The  tensile  data  for  the  Hercules  3501-6  epoxy  resin  discussed  by 
Browning  [32]  show  some  important  similarities  to  the  shear  data  shown 
in  Figure  7.  His  data  show  a  transition  temperature  above  which 
material  properties  are  degraded  greatly;  this  is  particularly  signifi¬ 
cant  at  high  moisture  contents.  This  transition  is  a  change  from  a 
"glassy"  solid  behavior  thigh  modulus) ,  to  a  "rubber"'1  material  behavior 
(low  modulus),  as  temperature  increases.  The  present  shear  data  (Figure 
7)  also  suggest  a  transition  temperature  over  the  entire  range  of 
moistures.  For  example,  the  zero  moisture  results  can  be  observed  as  a 
function  of  increasing  temperature.  At  100°C,  the  ultimate  octahedral 
shear  stress  has  decreased  as  has  the  ultimate  octahedral  shear  strain. 
But,  as  temperature  is  increased  further,  ultimate  strain  increases 
dramatically,  as  ultimate  stress  continues  to  decrease.  A  similar  re¬ 
sponse  is  observed  for  the  high  condition  as  temperature  increases. 
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Up  to  the  transition  temperature,  a  slight  decrease  in  ultimate 
strain  takes  place  as  the  ultimate  stress  drops  to  nearly  half  the  room 
temperature  value.  Temperature  has  a  softening  effect  (lowering  the 
modulus' ,  while  lowering  the  strength,  nut  the  material  still  has  a 
"giassv"  behavior,  as  seen  in  the  decreasing  ultimate  strain. 

'•.'hen  there  is  no  moisture  present,  low  temperature  ultimate 
stresses  and  ultimate  strains  from  the  shear  data  are  much  higher  than 
from  browning's  data  [32 1 .  At  the  high  temperature  of  160°C,  the  shear 
data  ultimate  stress  prediction  dips  below  Brown  Inc.  ’  -•  prediction,  but 
the  ultimate  strain  prediction  is  nearly  five  time;,  higher  than 
Brown ing ' s . 

Then  moisture  i  ;  added ,  the  trend  in  ultimate  j tresses  and  strains 
in  the  .-.hear  test  data  and  Browning’s  data  ire  similar .  Both  nets  of 
data  show  the  pi  a.  it Lc icing  effect  of  high  temperature  and  moisture. 

•vernl 1 ,  the  present  solid  rod  shear  data  show  much  more  pronounced 
inelastic  behavior  it  all  ootid  i  t  ions  than  the  t  ana  i  1 e  Tara  of  Browning. 
This  leads  to  higher  ultim.it. ■  stress  and  ttr.iin  prodi  it  ions. 

’-.'hen  comparing  the  oct  lhedral  shear  curves  <'i  '-ir.ure  b,  generated 
from  tensile  test;  data  i  }i;  the  present  invest  Lgat  nr,  to  Browning’s 

r,  'suits,  r 'lulus  ’-a  Ilf.  '  are  -oett  te  he  enis/  sLightlv  higher.  Again  at 
'.lgt  t  emn*.  r  1 1  •  ir  •  ■  inti  moisture  c  md  i  t  ions ,  the  present  tensile  data  snow 

i 

I  ! 

Lower  nodui  i  vaiue  .  ami  ltiwer  ultimate  stress  anti  strain  vacees.  The 
presence  '  a  gl  in  transition  is  on  1  v  hinted  at  in  tiie  present  data, 

•' v  r.he  ilii.ht  incre.i,.e  in  :  i  t  i  cut  t  e  strain  when  temperature  is  increased 
from  room  temperat  ure  to  !  i’  at  the  ..lturated  moisture  tontht  ion. 
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Overall,  the  highly  elastic  and  low  ultimate  strain  behavior  of  the 
tensile  data  follows  Browning's  tensile  data  more  closely  than  the 
shear  data  do.  This  is  as  ecpected  since  the  data  resulted  from 
similar  tensile  test  methods. 

Based  upon  the  above  comparisons,  the  shear  data  of  Figure  7 
appear  to  be  reasonable,  but  they  are  enough  unlike  the  tensile  data 
of  Figure  6  and  Browning's  results  [32]  to  require  more  experimental 
evidence.  More  precise  and  comprehensive  shear  data  tor  temperature 
and  moisture  variations  need  to  be  generated,  while  considering  the 
molecular  structure  and  chemistry  involved  in  the  resin  system. 

Fiber  properties  were  mainly  obtained  from  Hercules  [33]  and 
Owens-Corning  [ 34 ]  literature,  for  the  AS-graphite  fiber  and  the  52- 
glass  fiber,  respectively.  Shear  moduli  and  shear  stress  allowables 
necessary  for  stiffness  calculations  and  the  fiber  failure  criterion 
were  easily  found  for  the  S2-glass  (assumed  isotropic).  For  the  graph¬ 
ite  fiber,  in-plane  and  longitudinal  shear  tests  required  to  determine 
these  properties  are  not  commonly  performed  on  a  fiber,  due  to  its 
small  diameter.  Instead,  a  shear  modulus  value  was  calculated  in  the 
transverse  plane  of  symmetry  from  the  respective  values  of  Young's 
modulus  and  Poisson's  ratio.  It  was  necessary  to  estimate  values  of 
longitudinal  shear  modulus  and  longitudinal  shear  ultimate  stress  from 
data  for  other  similar  graphite  fibers  [35].  The  values  were  5.0  Msi 
and  225  ksi  respectively,  as  shown  in  Table  3.  The  transverse  shear 
stress  being  relatively  unimportant  because  of  the  inability  of  the 
analysis  to  directly  apply  transverse  shear  loads,  a  value  of  25  ksi 
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transverse  shear  allowable  was  chosen.  Microscopically  the  graphite 
fiber  is  much  Like  a  composite  due  to  the  longitudinally  arranged 
graphite  crystal  structure.  Thus,  a  low  transverse  allowable  shear 
stress  is  a  viable  approximation. 

Properties  of  the  constituent  materials  are  presented  xn  Table  3, 
with  estimated  values  denoted  bv  an  asterisk  (*).  Two  sets  of  epoxy 
matrix  material  properties  at  room  temperature,  dry  conditions  are 
also  shown,  generated  from  tensile  and  shear  experimental  data  (Figures 
6  and  7 ) . 


ticuutw 


SECTION  5 


NUMERICAL  RESULTS 

Test  cases  analyzed  using  the  longitudinal  shear  loading  version 
of  the  generalized  plan  strain  micromechanics  program  indicate  that 
the  new  analysis  performs  all  of  the  same  operations  as  the  old 
version,  with  the  addition  of  a  shear  loading  capability.  As  theory 
predicts,  there  is  no  coupling  between  shear  and  normal  loads  during 
elastic  load  increments,  but  during  increments  beyond  the  elastic  limit 
the  Prandtl-Reuss  flow  rule  is  in  effect  and  coupling  does  occur.  When 
the  deviatoric  stress  tensor  is  non-zero  (during  inelastic  increments), 
the  Prandtl-Reuss  flow  rule  causes  each  elemental  stiffness  matrix  to 
become  fully  populated.  This  in  turn  forces  each  stress  term  to  be 
dependent  on  each  strain  term.  Thus,  a  shear  stress  can  then  induce  a 
small  normal  strain,  which  is  impossible  in  elasticity  theory,  but 
entirely  possible  in  plasticity  theory.  In  fact,  all  stress  components 
have  an  effect  on  all  strain  components  if  the  deviatoric  stress 
tensor  is  fully  populated.  Looking  back  to  Eq .  (6),  it  can  be  seen 
that  the  inelastic  constitutive  equation  governs  the  coupling. 

The  additional  strain  induced  during  inelastic  increments  is 
due  to  the  two  deviatoric  stress  terms  within  each  term  of  that 
constitutive  equation.  This  effect  is  small,  but  significant  enough 
to  warrant  its  inclusion.  This  inelastic  strain  contribution  has  been 
observed  to  be  up  to  6  percent  of  the  total  strain.  To  demonstrate 
that  the  Longitudinal  version  of  the  micromechanics  analvsis  provides 
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accurate  and  useful  results,  two  types  of  examples  have  been  devised. 
The  first  uses  the  constituent  material  properties  to  predict  results 
using  longitudinal  shear  loading  analysis  which  are  then  compared  to 
actual  solid  rod  torsional  shear  data,  for  two  different  composite 
materials.  The  second  example  predicts  laminate  behavior  using  the 
longitudinal  shear  loading  analysis  in  conjunction  with  a  laminate 
point  stress  analysis. 

Comparisons  of  Analytical  Predictions  with  Solid  Rod  rsion  Test  Data 

Solid  rod  shear  data  had  been  generated  [27]  fo  ^mposite  spec¬ 

imens  as  well  as  for  the  previously  discussed  solid  rod  torsion  tests 
of  the  epoxy  matrix.  Fibers  used  in  the  composite  specimens  were 
the  Hercules  AS-granhite  fiber  [33]  and  the  Owens-C  -  -ning  "-glass 
fiber  [34].  Hercules  3501-6  epoxy  resin  was  used  as  the  matrix  system. 
These  composites  will  henceforth  be  referred  to  as  OR/EP  and  GL/EP, 
respectively,  for  the  graphite  and  glass  fiber  composites. 

Using  the  properties  evaluated  from  the  previously  discussed 
solid  rod  torsion  tests  of  the  matrix  system  alone  (Section  4), 
the  micromechanics  analvsis  predicted  the  composite  shear  stress- 
shear  strain  behavior.  Because  composite  tests  of  GR/EP  and  GL/EP 
were  performed  in  the  same  program  as  the  matrix  solid  red  torsion 
tests,  they  provided  the  logical  source  of  data  comparision.  These 
comparisions  are  shown  in  Figures  8  through  12  for  GR/EP,  and  Figures 
13  through  17  for  GL/EP.  Tests  were  performed  at  room  temperature, 
100°C,  160°C,  at  two  moisture  conditions,  i.e.,  dry  and  fully  saturated 
(6.75  percent  moisture  by  weight).  Although  data  were  also  available 
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Figure  10.  Shear  stress-shear  strain  curves  for  graphite/epoxy  at  100°C,  dry  conditions, 
comparing  micromechanics  predictions  (shown  in  triangles)  with  solid  rod 
torsion  test  data  (shown  in  solid  line) 
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Figure  16.  Shear  stress-shear  strain  curve 
by  weight)  conditions,  comparin 
with  solid  rod  torsion  test  dat 


55 


for  fully  saturated  composite  specimens  tested  at  160°C  for  both  the  GR/EP 
and  the  GL/EP,  the  analysis  showed  failure  initiation  occuring  before 
preconditioning  was  completed.  Thus,  these  cases  were  emitted.  This 
signifies  the  occurrence  of  microcracking  before  the  composite  is 
loaded  mechanically,  solely  due  to  thermal  and  moisture  loads.  Thus, 
not  only  do  temperature  and  moisture  increases  lower  the  matrix 
properties,  thev  induce  stresses  sufficiently  high  to  cause  failure. 

As  Figures  8  through  17  show,  in  most  cases  the  predicted  composite 
stress-strain  curves  were  in  close  agreement  with  experimental  results. 

At  present,  the  longitudinal  shear  micromeehunics  analysis  is 
only  operational  up  to  first  element  failure,  which  signifies  crack 
initiation  or  some  other  local  failure  on  the  micro  scale.  The 
present  lack  of  a  crack  propagation  capability  prevents  the  analysis 
from  predicting  the  actual,  failure  strength  of  a  composite;  but  it 
does  predict  inelastic  behavior  and  the  initiation  ->f  failure.  Thus, 
each  of  the  Figures  8  through  17  show  only  the  initial  portion  of  the 
complete  shear  stress-shear  strain  curve,  but  enough  to  reveal  the 
accuracy  of  the  prediction  technique. 

The  analysis  exhibited  excellent  agreement  with  experimental 
values  or  composite  modulus  for  the  GR/EP  specimens,  along  with  good 
predict  ion.'  of  viol d  strength.  This  was  a  significant  improvement 
over  the  results  obtained  when  using  tensile  data  to  generate  the 
ortuhedra  1  '.hear  stress-octahedral  shear  strain  curves  for  the  matrix 
materia  I  .  lite  tensile  lata  exhibited  much  less  nonlinearity,  as  indi¬ 
cated  in  Figure  n .  \1 though  the  use  of  data  generated  from  these 
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tensile  tests  caused  the  micromechanics  analysis  to  predict  low  values 
tor  first  element  failure  and  low  yield  predictions,  initial  modulus 
predictions  were  quite  accurate,  as  expected  since  the  initial  slopes 
in  Figures  6  and  7  are  similar.  Therefore,  the  tension  data  cannot  be 
refuted,  but  the  results  suggest  that  highly  inelastic  matrix  material 
behavior  does  occur  when  high  shear  loads  are  present.  Thus, 
subsequent  data  comparisions  in  this  section  wiil  employ  the  solid  rod 
torsion  test  data.  In  future  investigations,  the  physical  signigicance 
of  this  difference  between  the  two  test  methods  wiil  deserve  some 
attention. 

The  slight  but  consistent  overprediction  by  the  analysis  of  the 
experiment's  stress-strain  curves  for  GR/EP  will  be  noted ,  while  the 
results  for  GL/EP  are  more  scattered.  It  is  assumed  that  the  large 
modulus  difference  between  fiber  and  matrix  is  _o  blame  for  the  slight 
overpredictions.  Considering  Table  3  further,  the  contrast  between 
moduli  of  fibers  and  matrix  are  clearly  evident:  the  longitudinal  nor¬ 
mal  and  shear  moduli  of  both  fibers  are  much  higher  than  that  of  the 
epoxy.  For  future  consideration  it  will  also  be  noted  that  the  trans¬ 
verse  modulus  differences  from  fiber  to  matrix  are  much  higher  for 
GL/F.P  than  GR/EP.  A  general  statement  can  be  made  concerning  these 
differences:  the  higher  the  modulus  difference  between  fiber  and  ma¬ 

trix,  the  higher  the  stress  concentrations  in  a  composite  resulting 
from  an  applied  load  (loading  according  to  the  particular  direction  of 
those  moduli).  Large  modulus  difference  also  signifies  that  any  flaw 
in  the  composite  will  magnify  the  already  present  stress  concentrations , 
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thus  degrading  the  composite  properties.  The  analysis  idealizes  a 
composite  materia]  by  assuming  perfectly  homogeneous  materials  (neg¬ 
lecting  all  flaws),  and  perfect  bonding  between  fiber  and  matrix. 
Because  actual  composite  materials  have  ir.anv  flaws  in  the  form  of 
voids,  extraneous  inclusions,  and  debondings,  to  mention  the  most 
common,  the  present  analysis  is  expected  to  overpredict  slightly  in 
most  cases. 

Consider  the  drv  CR./EP  results  in  Figures  8,  10,  and  12.  Tire 
overpredic tion  error  decreases  as  temperature  increases;  at  160°C  the 
predicted  modulus  and  the  experimentally  measured  modulus  are  essen¬ 
tially  equal.  The  increased  temperature  causes  softening  of  the 
matrix  material,  thus  relieving  high  stress  concentrations.  Any  flaws 
and  imperfect  fiber-matrix  bonds  that,  occur  in  the  composite  affect  the 
total  behavior  of  the  composite  less,  thus  the  analysis  predicts  more 
accurately.  This  hypothesis  is  reinforced  by  a  studv  of  the  octahedral 
shear  stress  cent  sirs,  as  will  he  discussed  later. 

The  larger  errors  in  the  predictions  for  saturated  moisture  con¬ 
ditioned  0R/F.P  specimens  (.Figures  d  and  11)  can  be  attributed  to  the 
stress  intensity  caused  bv  high  moisture  dilation,  i'll  i  le  temperature 
softens  the  matrix  greatly,  moisture  causes  only  small  degradation  of 
material  stiffness,  is  shown  in  the  octahedral  shear  stress-octahedral 
shear  strain  curves  in  Figure  7.  The  higher  stresses  intensify  the 
detriment  il  effects  of  flaws,  causing  a  discrepancy  between  experiment 
and  theory.  Again,  the  trend  toward  lower  errors  .is  the  temperature  is 
increased  will  he  observed. 
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Although  trends  are  predicted  correctly  tor  the  Gl./EP  composites, 
there  is  a  significant  difference  between  the  anal  vs  is  and  the  actual 
experimental  data.  Figures  13,  Is  and  16  show  comparisons  at  room 
temperature,  dry,  room  temperature,  saturated  and  100°C,  saturated 
conditions,  respective! v .  These  examples  show  reasonably  close  agree¬ 
ment  between  theory  and  experiment.  Large  variations  are  observed  for 
both  of  the  elevated  temperature,  drv  test  conditions,  however,  as 
shown  in  Figures  L5  and  !7.  Because  the  matrix  at  elevated  temperature 
is  highly  nonlinear  and  very  low  in  strength,  the  prediction  is  made 
more  difficult  than  for  the  other  test  cases.  Also,  the  probability 
of  experimental  error  increased  at  the  highest  ter...  out  condition 
( 160°C)  due  to  the  limitations  in  temperature  c-v-.tr.,!  nailable  at  the 
time  of  the  testing,  as  discussed  in  Sect  ion  4.  ’  .  .•.mi  ding  the 

results  in  Figure  17  for  the  above  reasons ,  tin  <  1 :  :  or  EP 

represent  underpredict  ions .  As  temperature  in  t.o  under¬ 
predictions  grow  for  each,  separate  moisture  c.  mi:  •  .  :■<  cause  of 

tiiis  error  is  unknown  at  this  time;  a  stud'-  .  ■;  t:.«  o  ;  -ntours 

offered  no  explanation. 

Typical  behavior  of  a  shear  specimen  on  the  mi c :  •  scale  can  be 
seen  clearly  in  the  contour  plots  produced  t  rom  tin  mioromech.-mi.es 
results  of  dry  GR/EP  specimens  tested  at  room  temperature.  Octahedral 
•near  stresses  due  to  curing  ranged  up  to  75  percent  ot  the  elastic 
limit,  as  shown  in  Figure  18a.  Octahedral  shear  stress  is  an  indicator 
ot  distortional  energy  in  the  material ;  the  highest  values  are  seen  to 
occur  in  regions  of  close  fiber  spacing,  directly  between  the  fibers. 


a)  Octahedral  shear  stress 
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c)  Octahedral  shear  stress 
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stress,  T  =  4.5  ksi 
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Figure  18.  Contour  plots  of  octahedral  shear  (normalized  with  respect 
to  matrix  vield  stress,  4.18  ksi),  maximum  principal,  min¬ 
imum  principal,  interface  normal,  and  interface  shear  stress 
within  a  graphite/epoxv  solid  rod  torsion  specimen  at  room 
temperature,  dry  conditions. 


i)  Interface  shear  stress 
(ksi)  after  cooldown 
from  177°C  cure  temper¬ 
ature  to  21°C 


k)  Interface  shear  stress 
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j)  Interface  shear  stress 
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Figure  18.  (continued)  Contour  plots  of  octahedral  shear  (normalized 
with  respect  to  matrix  vield  stress,  4.13  ksi),  maximum 
principal,  minimum  principal,  interface  normal,  and  inter¬ 
face  shear  stress  within  a  graphite/epoxv  solid  rod  torsion 
specimen  at  room  temperature,  drv  conditions. 
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w:no  ••  e  1«  t.  \  if  .u  Figure  i  8h  )  and  eventual  iv  lirst  failure  (Figure 
lb,.  :  .iti,!  .  'it-  2. cm  contour  value  correspond  ing,  tj  i  irst  element 

•  .1 1  lur.-  r.'incH  tits  an  ootahe.iral  .shear  stress  of  !  1 .  xsi.  In  Figure 
~ i.ita  incic.ite  vielu  at  --3  ksi,  although  the  octahedral  shear 

.  :h»  ur  :  lot  ■  :  ’ure  !  Si:  >  shows  near  1’  2'<  ent  of  the  matrix 

•  *ave  •. roadv  viold*..:,  revc-ai  ing  the  :  act  that  local  matrix  vielding 
■  >  :  •  tic  are  the  ;  rap  os  i  1 1*  as  a  whole  begins  to  show  any  sign  of 

•ton  linear  response.  Fven  the  predicted  shear  s  tru-ss -'.hoar  str  tin 
resp.--.se  ’.shewn  bv  tr.e  triangles  in  Figure  8)  remains  nearlv  linear 
a •  • :  1  past  the  point  '  it  sr.  .  lemon c  vie’ ding. 

Vhile  the  ec  •  ahedra  i  ..  ear  stress  failure  criterion  predicts 
failure  direct.'  •  between  fibers,  a  c.ixinun  stress  criterion  would  dis¬ 
agree.;  ;u  ;h  tensile  stresses  ro  seen  to  occur  near  the  fiber  in 
Figure  i(5e.  High  .-m  r  ■:  •.  .■  ,  tresses  are  seen  to  .«o.  ur  at  the 

interface  along  the  ••.  :  '.xretrv  ■.  Figure  is:  •  .it  the  fiber-matrix 

interface,  indicating  *  p  •:  er.t  .  al  failure  ot  t!.<  inter:  i.-c  bond .  The 
octaher.ri'  drear  mru  ...  r: :  Ion  indicates  i  pr  c  !•  i  ntor:  iber 
r..i  cr.it  rack  'all  ire,  .mi.,:,  com  id  propagate  to  tin  t .  rface,  or  merely 
cause  a  split  direct  Ip’  between  adjacent  f  ibei  w ,  h.  t  in  criteria  predict 
failure  in  the  same  general  region  or  close  fiber  spacing,  but  when 
dealing  with  an  actual  composite,  the  actual  failure  mode  will  also  be 
influenced  by  the  particular  flaws  in  that  region.  For  example,  a 
iec.il  region  u  poor  '  iber/m.atrix  bonding  would  be  opened  by  the  high 
t-nsiln  stress  shown  in  Figure  I8e . 
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Cooldown  from  the  cure  temperature  induces  high  compressive 
normal  stresses  at  the  interface,  due  to  the  contraction  of  matrix 
material  around  the  fiber  (Figure  18g) .  The  contraction  of  the  matrix 
material  in  conjunction  with  the  constant  displacement  boundary  con¬ 
ditions  moves  the  boundaries  together  slightly,  enough  to  compress  the 
small  volume  of  material  near  the  axes  of  symmetry.  This  normal  stress 
is  unaffected  by  subsequently  applied  longitudinal  shear  loads,  and 
only  slightly  increased  by  contributions  due  to  inelastic  material 
response  at  higher  Levels  cf  the  applied  shear  stress,  as  shown  in 
Figure  18h.  The  reasons  for  this  are  that  there  is  no  coupling  between 
normai  and  shear  stresses  during  elastic  behavior,  and  that  when 
coupling  does  occur  (when  the  matrix  behavior  is  inelastic),  the 
deviatoric  stress  censor  is  so  sparsely  populated  that  the  normal 
stress  is  only  increased  by  about  three  percent. 

Due  to  the  additional  longitudinal  shear  stress  components,  the 
present  analysis  must  include  both  in-plane  and  out-->f-plane  shear 
stresses  in  the  interface  shear  stress  contour  plots,  in  contrast  to 
previous  treatments  (15,  30,  381.  Now  the  shear  stress  contour  repre¬ 
sents  only  the  magnitude;  the  direction  of  the  interface  shear  stress 
can  vary  from  position  to  position  along  the  fiber  in  any  direction  on 
the  fiber  surface. 

Interface  shear  stresses  developed  during  cooldown  (Figure  18i) 
are  only  25  percent  of  the  shearing  stress  at  failure  (Figure  181). 

It  will  be  noted  that  the  curing  stresses  become  purely  hydrostatic  on 
the  horizontal  and  vertical  axes  of  symmetry  and  at  45°,  due  to 


geometric  symmetry.  Viewing  Figures  1 8i ,  j,  k.  and  1  i-onsecut  ive  Lv 
shows  that  the  initial  curing  residual  shear  stress  small  relative 
to  the  load-induced  shear  stress,  which  eventually  initiates  failure. 

The  GR/EP  under  a  test  condition  of  room  temperature,  when  mois¬ 
ture-saturated  (6.75  percent  moisture  by  weight),  reveals  very  high 
octahedral  shear  stresses  and  interface  normal  stresses  due  to  mois¬ 
ture  dilatation,  as  seen  in  Figures  19a  and  19b.  Although  thermal  con¬ 
traction  during  cooldown  counteracts  the  subsequent  moisture  dilatation, 
the  extent  of  moisture  dilatation  far  overshadows  the  thermal  contrac¬ 
tion.  The  majority  or  matrix  material  is  already  inelastic  before 
loads  are  applied,  (Figure  L9a)  which  reduces  the  subsequent  load- 
carrying  capability,  while  the  moisture  also  softens  the  matrix, 
thus  reducing  its  ultimate  strength.  Figures  19c  and  1 9d  show  the 
increased  inelastic  behavior  that  entirely  envelopes  the  matrix  at 
failure.  The  moisture  loading  governs  the  contours,  i.e.,  they  remain 
more  nearly  symmetric  than  in  the  Loading  condition  without  moisture 
(Figure  I8d) .  Again  a  high  ;  yerface  normal  stress  is  present 
i  Figure  19e),  while  a  '  a  ’■efface  shear  stress  (Figure  19f)  con¬ 
tributes  to  the  high  octahedral  shear  stress  state  shown  in  Figure  19d . 
i he  maximum  and  minimum  principal  stress  contours  (Figures  19g  and  1 9h ) 
have  opposite  results  than  in  the  room  temperature,  dry  specimens,  due 
to  the  opposite  effects  of  moisture  and  curing  stresses.  The  minimum 
principal  stress  shows  a  highest  absolute  value  at  45°  from  the  x-axis 
of  symmetrv  at  the  interface.  Because  the  absolute  value  of  the 
maximum  principal  stress  is  higher  than  the  minimum  principal  stress. 


failure  would  be  predicted  to  occur  at  the  fiber-matrix  interface  near¬ 
est  the  x-axis  if  the  maximum  stress  criterion  were  used.  A  different 
result  is  again  shown  by  the  octahedral  shear  stress  criterion  (Figure 
I9d),  which  predicts  an  interfiber  failure,  although  it  is  near  the 
same  area. 

The  elevated  temperature  ( L00°C)  GR/EP  specimen  under  dry  condi¬ 
tions  shows  a  relatively  less  intense  octahedral  shear  stress  state 
due  to  curing  (Figure  20a)  than  the  room  temperature,  dry  specimen. 

This  occurs  not  only  because  of  the  smaller  temperature  excursion, 
but  also  because  of  the  lower  stiffness  properties  at  the  elevated 
temperature.  As  noted  earlier,  this  results  in  lower  stress  concen¬ 
trations  when  flaws  are  present  in  real  specimens.  Thus  the  assump¬ 
tion  o ;  no  flaws  and  perfect  bonds  becomes  less  severe,  and  the  anal¬ 
ysis  predicts  more  accurately.  Interface  normal  stresses  are  negli¬ 
gible  compared  to  the  magnitude  of  the  octahedral  shear  stresses 
occurring  at  failure  ^Figure  20b),  while  maximum  and  minimum  principal 
stress  contours  predict  maximum  absolute  values  at  the  same  point  as 
the  octahedral  shear  stress  contours:  between  the  fibers  on  the  x-axis 
of  syrametrv  (Figures  20c,  d).  Therefore,  an  elevated  temperature 
GR/F.P  specimen  would  tend  to  fail  within  the  matrix,  where  a  crack 
would  inflate  and  subsequent lv  begin  to  propagate. 

Because  the  GI./KP  at  room  temperature  and  6.75  percent  moisture 
has  a  smaller  fiber  volume  (50.5  percent)  than  the  GR/EP  at  similar 
conditions,  ind  because  the  fiber  properties  are  different  than 
granhite,  sLigntlv  different  results  can  be  expected  from  the 
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microana Lysis.  But  after  cooldown  from  the  raring  temperature  and 
moisture  condition  in*;,  eR/EP  and  u./EP  snow  very  similar  octahedral 
shear  stress  contours  isee  'igures  19a  and  21a).  Coefficients  of 
tuermal  expansion  varv  markedly  between  the  glass  fiber  and  the  graph¬ 
ite  fiber  (see  Table  ')■  The  graphite  fiber  expands  axially  when 
temperature  is  lowered,  and  contracts  in  the  transver-e  direction.  The 
glass  fiber,  being  isotropic,  contracts  in  all  material  directions 
when  temperature  is  1  'wered.  However,  the  matrix  thermal  contraction 
is  much  higher  than  char  of  either  fiber.  These  large  differences 
between  coeff  ;c lent  s  if  tuermal  expansion  from  fiber  to  matrix  cause 
stresses  to  :e  induced.  These  Jifferences  are  aighest  in  the  axial 
direction  of  tne  d\  and,  i  west  in  the  transverse  direction  of 

i.R/EP,  whi  ,  e  t  he  C'./T.P  is  in  between,  although  quite  high  (ree  Table  3) 
Cmsideritu;  this  and  the  fjc:  teat  :m  te  •i'*  r  is  present  in  the  Cl/EP 
bec.n-.se  of  d.  lower  :  irter  volume,  'iie  overall  thermal  and  moisture 
o  f  r  e  r  ••  ten;  :  ■  ne  equal  for  •'„ti  composites .  Therefore,  similar 
contours  a;e  ;ee..  for  both  nR  EP  and  CT./El'  in  figures  19a  and  21a. 

Because  o .  end  constraints  appl ied  in  thi  s  ana l vs is ,  hvgrothermal 
load.,  will  net  induce  lonp.J  t  i  n.-d  shear  stresses.  Considering  the 
t-.msv«rse  plane  .  >n  1  >,  tiie  overwiu- Imi  ng  moisture  dilatation  tends  to 
push  ini  ware  on  the  boundaries.  Since  the  boundaries  a  re  constrained 
to  nave  a  constant  displacement,  tiie  irea  of  highest  matrix  volume 
toward  tec  upper  right  corner  of  the  quadrant  will  be  slight Iv  in 
comnress ion .  while  the  area  of  small  matrix  volume  will  he  in  large 
tension,  as  s nowu  in  figure  21b.  Because  tiie  fiber  is  not  affected  bv 


.  k 


figure  21.  Contour  plots  of  octahedral  shear  (normalised  with  respect 
to  matrix  vield  stress,  2.88  ksi),  interface  normal  and  in¬ 
terface  shear  stress  within  a  glass/epoxv  solid  rod  torsioi 
specimen  at  room  temperature,  saturated  conditions  (6./J.M 
bv  we ight  '  . 
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this  moisture,  it  tends  to  hold  the  matrix  in  place  while  the  matrix 
tends  to  dilate  out  from  the  center.  Thus,  shear  stress  is  induced 
along  the  interface.  Because  matrix  dilatation  occurs  in  all  direc¬ 
tions,  and  because  of  geometric  symmetry,  only  a  hydrostatic  stress  is 
present  at  the  center  of  the  fiber/matrix  interface,  as  shown  in 
Figure  21c.  Constant  displacement  boundary  conditions  at  the  edge  of 
the  quadrant  plot  restrict  shear  strains  there,  and  again  only  a  hydro¬ 
static  stress  is  present.  This  is  shown  in  Figure  21c  by  the  absence 
of  shear  stress. 

As  the  applied  shear  stress  is  increased  until  the  first  element 
failure  occurs,  the  contours  propagate  in  a  manner  similar  to  previous 
cases.  But  comparing  Figures  2 l d  and  19d,  the  average  stress  in  the 
matrix  for  GL/EP  is  higher  iu.l  ksi)  than  that  for  CR/EP  (8.1  ksi) , 

The  matrix  in  GL/EP  has  yielded  to  a  greater  extent  than  in  GR/EP, 
signifying  more  iamage  has  taken  place,  and  more  energy  absorbed. 

Thus,  lower  fiber  volume  makes  a  more  uniform  stress  distribution  pos¬ 
sible  in  the  matrix,  i.e.,  the  stress  concentration  due  to  the  fiber 
decreases . 

The  brief  results  presented  here  show  good  correlation  with 
experimental  data,  suggesting  that  the  micromechanics  analysis  is 
valid.  These  results  cannot  be  considered  conclusive,  of  course.  A 
much  more  comprehensive  study  should  be  undertaken  in  the  future,  to 
reinforce  the  present  results. 
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Laminate  Analysis 

If  laminate  behavior  can  be  predicted  by  the  mieroanalysis  used 
in  conjunction  with  a  laminate  analysis,  further  value  and  accuracy 
of  the  micromechanics  analysis  can  be  demonstrated.  Longitudinal  shear 
stresses  are  indeed  a  necessary  consideration  in  laminate  analysis, 
and  are  capable  of  being  considered  in  the  present  micromechanics 
analysis.  The  use  of  micromechanics  analyses  was  restricted  to  uni¬ 
directional  laminates  in  prior  works,  when  no  longitudinal  shear  con¬ 
siderations  were  taken  into  account.  The  laminate  analysis  used  must 
be  compatible  with  the  restrictions  of  the  micromechanics  analysis, 
and  certain  assumptions  must  be  made.  Because  application  of  in¬ 
plane  shear  stress  (txv)  is  not  a  capabilitv  of  the  present  micro¬ 
mechanics  analysis,  a  two-dimensional  classical  laminated  plate 
point  stress  analysis  was  chosen  over  a  three-dimensional  finite 
element  laminate  analysis.  The  AC-3  laminate  point  stress  analvsis 
computer  program  is  operational  at  the  University  of  Wyoming  [ij. 
Basically,  it  describes  elastic  stresses  and  strains  induced  in  each 
ply  due  to  loads  applied  to  the  laminate.  These  ran  be  temperature 
and  moisture  Loads  as  well  as  mechanical  loads. 

Compatibility  of  the  two  analyses  and  the  scheme  of  arriving  at  an 
inelastic  stress-strain  curve  for  a  particular  laminate  is  described 
by  the  following  four  steps: 

1.  Calculate  ply  properties  from  the  mi cronechanics  analvsis. 

2.  Applv  the  laminate  point  stress  analvsis  to  the  prescribed 

lamirite,  loaded  in  the  prescribed  manner,  to  obtain  the  stress 


73 


state  induced  in  each  ply. 

3.  Use  the  micromechanics  analysis  to  analyze  each  unique  ply 

of  the  laminate,  holding  loads  in  the  same  ratios  as  the  laminate 
analysis  predicted.  This  will  result  in  inelastic  stresses  and 
strains  in  each  ply  when  loaded  to  high  levels. 

4.  Transform  the  stress  and  strain  states  back  to  laminate  coor¬ 
dinates  to  obtain  a  laminate  stress-strain  curve.  This  curve  can 
be  used  for  comparisions  with  laminate  experimental  data. 
Experimental  data  were  available  for  a  57.5  percent  fiber  volume, 

f ± 4 5  ]  .  GR/EP  laminate  at  four  combinations  of  temperature  and  mois- 
4s 

ture  conditions  [30),  viz,  room  temperature,  dry  (RTD) ;  room  tempera¬ 
ture,  one  percent  moisture  by  weight  (RTW) ;  elevated  temperature  (103 
°C),  dry  (etd) ;  and  elevated  temperature  y 1 0 3 °C),  one  percent  moisture 
'ETW).  The  conditions  will  hence  be  referred  to  using  the  abbrevia¬ 
tions  in  parentheses;  i.e.,  RTD,  RTW,  ETD,  ETW,  similar  to  the  notation 
of  Reference  [30].  Summarizing,  Reference  [30]  was  an  investigation  of 
compression  fatigue  properties  of  composites.  Static  compression  and 
compression  fatigue  tests  were  performed  on  materials  at  various  hy- 
grothermal  conditions.  Theoretical  predictions  of  failure  mechanisms 
were  studied  through  a  micromechanics  analysis  [1].  The  static  com¬ 
pression  experimental  data  for  the  [145]^s  laminate  were  chosen  for  the 
purpose  of  comparing  them  with  theoretical  predictions.  Both  constit¬ 
uent  materials  were  the  same  as  the  graphite  and  epoxy  proviously  de¬ 
scribed  in  Section  3. 


To  begin  the  analysis,  unidirectional  plv  properties  are 
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calculated  using  the  micromechanics  analysis,  to  be  consistent  with  the 
properties  utilized  later  when  inelastic  strains  are  to  be  calculated. 
Constituent  properties  for  the  AS-graphite  fiber  and  the  epoxy  matrix 
are  the  same  as  those  shown  in  Section  3.  Composite  properties  are 
shown  in  Table  4;  the  stiffness  values  correspond  closely  to  properties 
found  by  Northrop  Corporation  and  the  University  of  Wyoming  [30].  How¬ 
ever,  thermal  and  moisture  coefficients  of  expansion  differ  consider¬ 
ably.  It  is  important  to  be  mutually  consistent  in  the  present  devel¬ 
opment  from  one  analysis  to  the  other,  i.e.,  properties  used  initially 
in  the  laminate  analysis  must  correspond  to  values  (see  Table  4)  used 
later  in  the  micromechanical  analysis.  Overall,  temperature  and  mois¬ 
ture  lower  all  the  properties  recorded.  Also,  note  the  small  value  for 
longitudinal  thermal  expansion  coefficient  at  the  elevated  temperature, 
due  to  the  negative  coefficient  of  thermal  expansion  of  the  graphite 
constituent. 

The  ply  properties  as  calculated  using  the  micromechanics  analysis 
are  used  in  the  laminate  point  stress  analysis,  where  ply  stress 
states  due  to  temperature,  moisture  and/or  applied  loads  are  obtained. 
In  applying  a  thermal  or  moisture  load  to  the  laminate,  the  laminate 
analysis  assumes  the  input  elastic  properties  to  remain  constant 
throughout  the  temperature  or  moisture  change.  Because  the  matrix  mod¬ 
ulus  is  actually  not  constant,  but  decreases  significantly  with  in¬ 
creasing  temperatures,  the  predicted  curing  stresses  will  be  higher 
than  they  actually  are.  To  more  accurately  estimate  the  curing  stres¬ 


ses,  an  effective  temperature  change  can  be  used  which  is  smaller  than 
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TABLE  4 


PROPERTIES  CALCULATED 
FOR  USE  IN 


FROM  MICROMECHANICS  ANALYSIS 
LAMINATE  ANALYSIS 


Predicted  Elastic  Properties  at 
Environmental  Conditions  Indicated 


Properties 

at 

Predicted 

First 

Failure 


PROPERTY 

RTD 

RTW 

ETD 

F.TW 

RTD 

EL  (Msi) 

18.700 

18.690 

18.563 

18.554 

18.409 

Et  (Msi) 

1.467 

1.450 

1.172 

1.147 

1.559 

GLT(Msi) 

0.918 

0.897 

0.611 

0.589 

0.441 

VLT  , 

0.2573 

0.2571 

0.2544 

0.2542 

0.2805 

a,  (  x  10_6/°C) 

0.477 

0.455 

0.166 

0.144 

0.531 

aT 

32.40 

32.35 

31.70 

31.70 

33.86 

dL  (  x  10  7%M) 

0.042 

0.041 

0.026 

0.025 

0.046 

bT 

1.082 

1.078 

1.031 

1.027 

1.184 

the  actual  temperature  change.  For  curing  from  177°C  to  21°C,  an 
effective  temperature  change  of  111°C  was  used  rather  than  the  actual 
156°C .  In  curing  from  177°C  to  103°C,  only  a  42°C  effective  temper¬ 
ature  change  was  used.  These  assumed  effective  temperature  change 
values  were  taken  from  Reference  [ 30] .  Actual  loads  to  be  applied  to 


model  stresses  induced  in  the  individual  plies  during  each  case  of 


curing  and  conditioning  are  shown  in  Table  5,  along  with  actual  load 
ratio  increments  due  to  the  subsequently  applied  axial  compressive  load. 
To  arrive  at  these  values,  the  laminate  analysis  uses  the  composite  ply 
properties  that  were  calculated  by  the  micromechanics  analysis.  The 
resulting  stress  state  is  calculated  in  both  ply  coordinates  and  lamin¬ 
ate  coordinates  by  the  laminate  analysis,  thus  eliminating  the  neces¬ 
sity  of  transforming  the  stress  state  by  hand  for  subsequent  input  to 
the  micromechanics  analysis. 
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TABLE  5 

LOAD  RATIOS  CALCULATED  FROM  LAMINATE 
ANALYSIS  FOR  USE  IN  MICROMECHANTCS 
ANALYSIS  FOR  [+45]  LAMINATE 

l  J  4g 


HYDROTHERMAL  LOADS 


STRESSES  (ksi) 

RTD  i-lll°C) 

RTW  (1%  M,-UI°C) 

ETD  f-42°C) 

ETW  (  1,1  M,-4  2°C) 

5, 

-4.65 

-3.23 

-1.41 

-().  i'i 

fS 

4 .  o  5 

3.23 

1.41 

0.33 

~xz 

0.0 

0.0 

0.0 

0.0 

APPLIED  LOAD  RATIOS  FOR  -1000  psl  AXLAL  COMPRESSIVE  STRESS 


STRESSES 

psi) 

RTD  I 

_ 

RTW 

ETD 

ETW 

1 

-912  I 

-912 

!  t 

-92S  ! 

-930 

X 

-  88 

-  38 

r  i 
r- 

1 

-  70 

L-  ■  X*  .  ■■ 

300 

500 

300  ! 

500 

The  micromechanics  analysis  is  now  assumed  to  model  a  single  uni¬ 
directional  ply  of  the  laminate  in  ply  coordinates  by  modeling  one 
quadrant  of  the  repeating  unit  cell,  as  depicted  in  Figure  22.  The 
assumption  that  boundaries  displace  uniformly,  was  described  earlier  in 
Section  3.  The  laminate  analysis  intrinsically  makes  this  assumption, 
since  classical  point  stress  analyses  do  not  include  interlaminar 
shear  stresses.  In  an  actual  composite,  certain  interlaminar  shear 
stresses  cause  shear  deformation  in  the  plane  transverse  to  the  fiber 
direction.  This  is  in-plane  deformation  for  the  raicromechanics  model, 
which  is  inadmissible  in  this  analvsis  due  to  the  constant  displacement 
boundary  conditions.  Therefore,  interlaminar  shear  stresses  ere  as¬ 
sumed  to  be  negligible  in  the  present  example.  Another  assumption  is 


jue  to  inelastic  ma- 


that  the  load  ratios  will  lunge  only  negligibly  c 
teriai  behavior.  One  way  to  check  this  assuror  ion  is  by  calculating 
the  composite  material  properties  c a  ply  while  it  is  loaded  into  the 
inelastic  range.  These  revised  material  properties  (see  last  column  of 
Table  4)  were  found  to  v.irv  onlv  slightly,  except  for  the  shear  modulus 
which  decreased  to  a  value  of  about  half.  >:■_  to  the  15°  ply  orienta¬ 
tion  of  the  laminate,  • he  major  loading  stress  is  longitudinal  shear 
loading.  At  this  high  inelastic  range  of  loading,  the  slope  of  the 
shear  stress-shear  strain  curve  is  much  less  than  the  elastic  shear 
modulus,  explaining  t he  low  composite  shear  modulus.  These  values  can 
he  used  in  tie  laminace  anal vs is  to  *'ind  new  road  ratios,  which  when 
compared  to  the  original,  load  ratios  show  ,.nly  small  errors  (Tatle  6). 
The  new  load  ratio.-  '•••or-'  calculated  f <i.v  :  taken  after  the  plv  was 


•  ids  N  .VXAi..  .-  1  '  r  L.l’.u'l  H.\  1  .0  .  ;C\  1  ’\  1,1.1 

ENT i KM  LOAD  EXCURSION 

1001 »  psi  Applied  .‘•::ia  L  Compressive  Stress 
STRESSES  (ps  ;;l  ~XOPERTIES  j  IllEl^\ST"~"i^FRTTES |  PERCENT  ERROR 

I  :  ~'h% 


11. 87 
4.10 
0.0 


axially  loaac-d  to  -JO  ksi  to  assure  maximum  accumulated  error  in  load 
ratios.  Notice  in  Taole  o  the  high  axial  load  and  longitudinal  shear 
load  while  the  transverse  normal  load  is  small.  The  9.1  percent  error 
accumulated  in  the  small  transverse  load  is  negligible  compared  to  the 
large  ixl.i!  normal  load  accumulation  of  less  than  one  percent  error. 

'he  '.cr.gitui:iu.il  snear  load  experienced  no  change,  which  is  expected 
in  tins  where  botii  load  direction  and  fiber  direction  occur  in  a 

manner  such  tii.it  tne  material  can  be  considered  specially  orthotropic. 
Because  there  is  no  coupling  between  normal  and  shear  effects  in  the 
special.lv  irti’.ci  rep :  v  case,  no  change  should  occur  in  shear  stress, 
even,  m  ciu  inelastic  region.  Therefore,  it  is  possible  to  conclude 
that  no  significant  error  in  theoretical  predictions  can  be  attributed 
to  variation  in  load  ratios  due  to  inelastic  behavior. 

With  the  above  assumptions  in  effect,  the  entire  curing  and  con¬ 
ditioning  history  of  each  pry  is  approximated,  as  well  as  the  effect  of 
the  applied  load  increments.  The  resulting  inelastic,  strains  at  each 
incremental  applied  load  are  in  ply  coordinates  and  must  be  transformed 
back  to  laminate  coordinates.  This  done,  longitudinal  stress  and  strain 
at  various  increments  for  eacn  condition  are  the  final,  results  needed  to 
compare  with  experimental  data.  Figures  23  through  26  show  the  predic¬ 
ted  incremental  stress-strain  points  compared  to  the  experimental  data 
curves.  The  complete  experimental  data  curves  to  failure  are  not  shown 
here,  although  thev  are  available  in  Reference  f  30]  . 

A  verv  slight  overprediction  error  is  seen  for  both  of  the  room 
temperature  conditions,  while  the  elevated  temperature  condition 
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Figure  24.  Comparison  of  room  temperature,  wet  (RTW)  experimental  data  [30]  for  a  ['45 
laminate  with  combined  micromechanical  and  laminate  analysis  predictions. 
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Figure  26.  Comparison  of  elevated  temperature,  wet  (ETW)  experimental  data  [30]  for  a 

[±45]48  laminate  with  combined  micromechanical  and  laminate  analysis  predictions. 
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predictions  are  in  excellent  agreement  with  experiment.  Although  the 
differences  between  theory  and  experiment  are  perhaps  not  significant, 
the  most  likely  explanation  for  slight  differences  is  based  on  flaws  in 
the  matrix  or  flaws  in  the  fiber-matrix  bond,  as  hypothesized  previously 
for  the  case  of  pure  shear  loading.  Again,  the  stress  concentrations 
caused  by  these  flaws  reduce  the  strength  of  the  material.  The  RTD 
test  case  (Figure  23)  does  show  this  discrepancy.  Any  softening  of 
the  matrix  should  reduce  the  stress  concentrating  effect  of  flaws. 

When  moisture  is  added  in  the  RTW  case,  (Figure  24)  it  serves  as  a  plas¬ 
ticizer  that  tends  to  negate  the  thermal  curing  stresses.  Now  the 
actual  composite  behavior  is  slightly  closer  to  the  predicted  curve  for 
a  composite  with  no  material  flaws.  This  theory  is  further  supported 
in  the  elevated  temperature  test  cases  (Figures  25  and  26),  where  the 
predicted  incremental  values  are  in  very  close  agreement  with  the  ex¬ 
perimental  data.  The  elevated  temperature  and  moisture  in  the  ETW  test 
case  (Figure  26)  again  softens  the  matrix,  reducing  the  effect  of  stress 
concentrators,  and  narrowing  the  gap  between  experiment  3nd  theory. 

Neglecting  interlaminar  shear  stresses,  calculating  load  ratios 
using  an  elastic  analysis,  and  assuming  an  effective  temperature  change 
were  other  possible  sources  of  error.  All  of  chese  effects  could  de¬ 
finitely  increase  the  stresses  in  the  matrix,  which  would  effectively 
increase  the  matrix  "damage"  consequently  lowering  the  overall  stiffness 
and  strength  of  the  actual  specimen  below  the  predicted  theoretical 


values . 
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The  reason  for  Che  excellent  accuracy  of  the  analysis  can  be  investi¬ 
gated  by  observing  the  loading  history  recorded  in  the  stress  contour 
plots.  Also,  insight  into  other  possible  sources  of  error  may  be 
attained.  Normalized  octahedral  shear  stress  contours  for  the  RTD  test 
case  are  shewn  in  Figure  27,  ranging  from  curing  stresses  to  the  failure 
stress  state.  Referring  back  to  the  original  step-by-step  description 
of  the  laminate  analysis  presented  earlier  in  this  subsection,  it  will 
be  noted  that  hygrothermal  effects  are  considered  in  separate  steps. 
First  there  are  stresses  induced  by  hygro thermally-influenced  consis- 
tuent  material  properties,  (see  Figure  27a).  Second,  other  plies  are 
influenced  by  hygrothermal  effects,  which  influence  the  stresses  in  the 
representative  ply  now  being  considered  (see  Figure  27b).  These  will 
be  termed  "ply  curing  effects"  and  "laminate  curing  effects",  respec¬ 
tively.  Because  Figure  27a  is  also  representative  of  stresses  in  a 
unidirectional  composite,  contrasting  Figures  27a  and  b  shows  curing 
stress  differences  between  a  unidirectional  laminate  and  a  [±45]  lam¬ 
inate.  But,  for  the  present  purposes,  it  is  important  to  note  the 
additional  curing  effects  of  the  ply  configuration  in  Figure  27b  and 
the  magnitude  of  them.  Any  material  defect  which  might  exist  in  a 
local  region  of  high  stress  will  cause  a  severe  stress  concentration, 
which  in  turn  can  cause  premature  matrix  yielding  and  failure.  This  is 
a  possibLe  reason  for  the  lower  stiffness  exhibited  by  the  experimental 
data  at  room  temperature  conditions.  Figure  27c  shows  the  amount  of 
matrix  that  has  already  yielded  locally  when  the  overall  composite  re¬ 
sponse  first  indicates  a  yield  point  experimentally.  While  the  tangent 
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Figure  27.  Contour  plots  of  octahedral  shear  (normalized  with  respect 
to  matrix  yield  stress,  4.18  ksi),  interface  normal  and  in' 
terface  shear  stress  within  a  graphite/epoxy  [ ±4 5 ] lam¬ 
inate  ply  at  room  temperature,  dry  conditions  (RTD). 
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Figure  27.  (Continued)  Contour  plots  of  octahedral  shear  (normalized 
with  respect  to  matrix  yield  stress,  4. IS  ksi),  interface 
normal  and  interface  shear  stress  within  a  graphite/epoxy 
[±451,  laminate  ply  at  room  temperature,  dry  conditions 


53 


88 


modulus  of  Che  yielded  matrix  is  still  close  to  the  initial  modulus, 
the  nonlinearity  is  not  apparent  in  the  composite  response.  Only  after 
the  yielded  region  has  propagated  upward  a  distance  equal  to  half  the 
width  of  the  quadrant,  do  the  experimental  data  begin  to  suggest  an 
elastic  limit.  When  the  micromechanics  analysis  predicts  first  element 
failure,  this  implies  the  initiation  of  a  crack.  To  prevent  this  dam¬ 
age  initiation,  the  composite  should  thus  not  be  loaded  past  the  stress 
represented  by  the  final  predicted  point  in  Figure  23.  A  more  compre¬ 
hensive  study  would  reveal  just  how  much  load  could  safely  be  added  be¬ 
yond  yield,  while  preventing  major  microcracking.  After  curing,  the  com¬ 
posite  experiences  large  normal  stresses  along  the  interface,  as  shown  in 
Figure  27e,  suggesting  a  high  probability  of  debonding  behavior.  There 
rre  also  fairly  high  shear  stresses  along  the  interface  in  this  as-cured 
condition,  which  grow  enormously  as  load  is  applied  and  increased  to 
predicted  iirst  failure  (Figure  27f  and  27h) .  This  shear  stress  is 
the  combined  effect  of  the  applied  longitudinal  shear  stress  and  the 
transverse  normal  stress,  which  are  acting  in  weak  directions  of  the 
material,  and  which  eventually  lead  to  failure.  It  will  also  be  noted 
in  Figure  27g  that  the  normal  interface  stress  has  been  reduced  by  the 
applied  loads.  This  is  further  evidence  of  how  applied  loads  can  re¬ 
lieve  some  of  the  stresses  induced  during  cooldown;  it  is  an  inelastic 
coupling  effect. 

The  ETW  test  case  involves  two  separate  forms  of  curing  and  con- 
titioning  stress.  First,  stresses  are  induced  independently  in  each 
ply  by  thermal  and  moisture  changes  due  to  the  difference  in  properties 


between  the  fiber  and  matrix.  Secondly,  since  the  individual  plies, 
exhibiting  the  aforementioned  hygrothermal  strains,  are  bonded  to  each 
other,  they  mutually  affect  each  other  because  of  their  differing  ply 
orientations.  The  combination  of  thermal  and  moisture  loads  has  far- 
reaching  consequences  because  of  the  opposing  effects  in  cooldown 
and  moisture  absorption.  The  temperature  drop  during  cooldown  causes 
the  matrix  to  contract  by  a  relatively  large  amount,  the  fiber  to  con¬ 
tract  transversely  by  a  lesser  amount,  and  the  fiber  to  actually  expand 
slightly  in  its  axial  direction,  (see  Table  3).  The  composite  axial 
strain  is  predicted  to  be  negative,  as  can  be  seen  by  the  positive  com¬ 
posite  thermal  coefficient  of  expansion  in  Table  4.  Figure  28a,  when 
compared  to  Figure  28b  and  28c,  indicates  that  the  octahedral  shear 
stresses  induced  by  curing  are  almost  fully  negated  by  the  absorption  of 
one  weight  percent  moisture  by  the  composite.  The  small  hygrothermal 
stresses  which  result  are  the  reason  for  the  low  average  laminate  hygro¬ 
thermal  stresses  shown  in  the  last  column  of  Table  5.  However,  the 
shifting  of  the  octahedral  shear  stress  contours  in  Figure  28c  due  to 
the  application  of  the  induced  loads  will  be  noted.  This  is  depicted 
more  clearly  at  first  failure,  as  represented  by  Figure  28d.  The  detri¬ 
mental  effect  of  unfavorable  temperature  and  moisture  combinations  is 
clearly  shown  in  comparing  the  octahedral  shear  stress  contours  at  first 
failure  of  the  ETW  test  case  compared  to  the  RTD  case.  An  octahedral 
first  element  failure  stress  of  approximately  11.6  ksi  (shown  normal¬ 
ized  in  Figure  27d  for  the  RTD  case  as  the  2.8  contour)  is  reduced  to 
5.2  ksi  in  the  ETW  case  (shown  normalized  in  Figure  28d  as  the  2.2 
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Figure  28.  (continued)  Contour  plots  of  octahedral  shear,  interface 
normal  and  interface  shear  stress  within  a  graphite/epoxy 
[±45],  laminate  ply  at  103°C,  1"  moisture  by  weight  (ETW) 
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contour).  This  is  a  significant  change,  considering  that  it  is  due  to 
a  relatively  modest  temperature  change. 

Interface  normal  stresses  induced  during  cooldown  and  one  percent 
moisture  absorption  are  relatively  small  tensile  stresses  (see  Figure 
28e),  but  subsequent  mechanical  loading  of  the  laminate  counteracts  the 
induced  tensile  stresses,  so  that  near  the  x-axis  of  symmetry  a  small 
negative  value  exists  at  first  element  failure  (Figure  2Sg) .  Small 
tensile  stresses  are  induced  near  the  v-axis  of  symmetry  at  first  ele¬ 
ment  failure,  but  are  small  enough  that  they  have  no  noticeable  effect 
on  subsequent  laminate  behavior.  Although  the  interface  shear  stress 
after  curing  and  moisture  conditioning  (Figure  28f)  is  small  (4  per¬ 
cent  of  the  ultimate  stress),  the  subsequently  applied  mechanical 
loads  result  in  shear  stresses  along  the  interface  which  are  very  high 
at  first  element  'allure  (Figure  28h) .  Although  the  interface  shear 
stress  is  high,  cue  octahedral  shear  stress  is  shown  to  be  highest  at 
the  outer  edge  of  the  model,  along  the  x-axis  of  symmetry’  (Figure  28d) . 
Failure  in  this  elevated  temperature,  one  percent  moisture  condition  would 
probably  occur  between  the  fibers  in  the  region  of  highest  octahedral 
shear  stress.  But  at  room  temperature,  dry  conditions,  higher  stresses 
would  be  induced  along  the  interface  (see  Figure  27d),  increasing  the 
probability  of  failure  of  the  fiber/matrix  bond.  A  resulting  hypothe¬ 
sis  is  that  because  this  failure  is  likely  to  occur  more  often  in  the 
RTD  laminate,  an  actual  laminate  at  these  conditions  has  a  higher 
probability  of  failure  initiation  due  to  inherent  flaws  in  the  mater¬ 
ial,  especially  along  the  fiber/matrix  interface.  Thus,  experimental 
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data  show  a  lower  initial  modulus  and  a  lower  overall  stiffness  than 
theory  predicts,  as  in  Figure  23. 


SECTION  6 


DISCUSSION 

The  establishment  of  a  generalized  plane  strain  finite  element 
formulation  that  includes  a  longitudinal  shear  loading  capability  has 
further  advanced  composite  materials  technology,  by  making  detailed  in¬ 
formation  concerning  the  micromechanical  shear  loading  response  of 
composite  materials  available.  Such  a  quasi-three-dimensional 

formulation  approaches  the  capability  of  a  true  three-dimensional 
analysis,  while  retaining  the  conciseness  of  a  two-dimensional  analysis. 
Although  the  present  generalized  plane  strain  analysis  only  approximates 
a  true  three-dimensional  analysis,  the  present  results  show  that  it  is 
quite  accurate,  and  therefore  valuable  in  many  potentia’  applications. 

The  capabilities  of  the  new  analysis  developed  in  the  present  study 
will  first  be  summarized.  Following  a  previous  formulation  [1],  an  ine¬ 
lastic  capability  is  realized  through  the  use  of  the  tangent  modulus 
approach,  with  revisions  due  to  the  addition  of  longitudinal  shear 
loading.  This  capability  enables  the  analysis  to  accommodate 
highly  inelastic  matrix  materials  such  as  annealed  aluminum.  The  finite 
element  formulation  was  adapted  by  making  extensive  revisions  of  the 
stiffness  formulation  and  Gaussian  elimination  solution  procedures.  A 
special  loading  technique  was  developed,  to  make  possible  a  decrease  in 
the  size  of  the  stiffness  matrix  storage,  increasing  the  efficiency  of 
the  program  considerably.  This,  combined  with  a  new,  more  efficient 
finite  element  model,  results  in  a  stiffness  matrix  which  has  the  least 
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possible  amount  of  wasted  space.  Even  with  this  improved  efficiency, 
however,  the  added  longitudinal  shear  loading  capability  quadruples  the 
size  of  the  stiffness  matrix,  leading  to  the  conclusion  that  if  a  certain 
problem  does  not  necess itate  longitudinal  shear  considerations,  the 
previous  formulation  [1]  which  does  not  include  a  longitudinal  shear 
capability,  should  be  employed.  Although  the  two  analyses  could  be 
combined  into  one,  and  a  selection  process  built  into  the  resulting  com¬ 
puter  program,  it  would  essentially  be  just  that,  i.e.,  two  separate 
programs  linked  together,  with  no  special  advantage  being  gained. 

Consistent  with  the  previous  analysis  [1],  another  capability  of 
the  new  version  is  combined  loadings;  that  is,  increments  of  tempera¬ 
ture,  moisture,  and  five  separate  applied  mechanical  stresses  can  be 
applied  in  any  order,  simultaneously,  or  in  any  combination.  The  in¬ 
clusion  of  temperature-  and  moisture-dependent  matrix  material  proper¬ 
ties  provides  an  added  dimension  in  modeling  real  physical  behavior. 

The  program  is  thus  a  highly  versatile  analytical  tool,  with  many 
applications  yet  to  be  explored,  suggesting  that  future  work  should 
follow  the  directions  discussed  below. 

Further  verification  of  this  new  analysis  is  not  requisite  or 
urgent,  the  existing  verifications  being  reasonably  conclusive,  but 
such  work  is  suggested  while  additional  developments  are  being  under¬ 
taken.  A  comprehensive  three-dimensional  study  of  micromechanics 
should  be  performed  to  establish  the  amounts  of  error  incurred  by  the 
special  assumptions  of  generalized  plane  strain.  Such  a  three-dimen¬ 
sional  program  is  now  operational  [36].  Usefulness  of  the  present 
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analysis  can  thus  be  shown,  as  well  as  those  cases  where  it  is  most 
accurate. 

With  the  longitudinal  shear  loading  capability,  interfacing  with  a 
laminate  point  stress  analysis  is  made  possible,  enabling  the  investiga¬ 
tion  of  complicated  laminates.  The  results  of  the  comparisons  between 
theory  and  experiment  for  the  [±45]^g  laminate  presented  here  are  en¬ 
couraging,  although  limited  to  one  simple  laminate.  However,  because 
this  particular  laminate  is  dominated  by  a  longitudinal  shear  loading 
mechanism,  it  is  assumed  that  further  investigation  of  various  other 
laminates  will  again  show  good  agreement  with  experiment. 

The  combined  micromechanics/laminate  analysis  predicts  only  a 
slightly  nonlinear  stress-strain  curve  up  to  the  point  where  first 
element  failure  is  predicted  to  occur.  In  the  present  formulation,  this 
is  the  termination  point  of  the  analysis.  If  a  crack  initiated  at  the 
firs',  element  failure  is  allowed  to  grow,  perhaps  an  entire  stress- 
strain  curve  could  be  modeled,  and  ultimate  stresses  predicted.  This 
crack  propagation  capability  has  been  developed  in  another  micromechan¬ 
ical  analysis  program  [37],  which  could  be  closely  followed  in  incorpor¬ 
ating  a  similar  capability  in  the  present  program.  This  would  provide 
a  capability  to  more  completely  model  the  behavior  of  laminates  over 
their  entire  loading  range. 

The  trend  toward  adhesive  bonding  of  composites  rather  than  using 
mechanical  fasteners  suggests  further  analysis  of  bonding  behavior. 
Longitudinal  shear  loading  and  viscoelastic  considerations  would  be 
desirable  in  a  study  of  bonds  in  conjunction  with  laminates.  Because 
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a  nonlinear  viscoelastic  capability  is  already  operational  in  an  exist¬ 
ing  raicromechanics  formulation  [38],  the  present  analysis  could  be 
revised  to  incorporate  this  time-dependent  response,  and  thus  handle 
special  cases  of  generally  orthotropic  laminates  bonded  together  with 
adhesives,  all  exhibiting  viscoelastic  response. 

Implications  of  the  assumption  of  a  perfect  fiber /matrix  bond  have 
been  recognized  by  previous  authors  [1].  This  is  now  especially  impor¬ 
tant  when  high  longitudinal  shear  loads  are  applied.  Results  presented 
here  for  shear  test  comparisons  and  laminate  comparisons  with  actual 
experimental  data  both  show  an  overprediction,  thought  to  be  caused,  at 
least  in  part,  by  this  assumption  of  a  perfect  fiber /matrix  interface 
bond.  The  differences  are  seen  to  decrease  rapidly  with  increasing 
temperature  and  moisture  conditions,  at  which  the  plasticized  matrix  is 
more  readily  able  to  accommodate  the  local  stress  disturbances  caused 
bv  the  debonding  in  the  actual  material.  Statistical  models  of  flaw 
distributions  may  also  be  developed  in  the  future,  so  that  material 
properties  can  be  adjusted  to  account  for  such  damaging  conditions. 

This  would  possibly  improve  the  already  good  predictions  of  the  present 
analysis . 

Such  further  developments  of  the  analysis  presented  here  should 
Lead  to  the  continued  advancement  of  the  design  of  high  performance 
composite  structures.  As  such  analyses  develop,  design  parameters  for 
composites  will  become  more  useful  and  well-known.  Composites  will  in 
turn  become  the  logical  material  choice  for  many  applications,  present¬ 
ing  an  important  alternative  to  the  age  old  use  of  homogeneous  metals. 
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APPENDIX  A 


ISOTROPIC  CONSTITUTIVE  EQUATION 
FOR  ELASTOPLASTIC  BEHAVIOR 


The  Prandtl-Reuss  flow  rule  is  employed  to  characterize  elasto- 
plastic  behavior  by  relating  the  plastic  strain  tensor  to  the  devi- 
atoric  stress  as  depicted  in  Eq .  (A-l) : 

eij(P)  =  'si j  (A_1> 


where  j  is  the  plastic  strain  tensor,  \  is  a  positive  scalar,  and 
s^j  is  the  deviator ic  stress  tensor.  From  this,  a  constitutive 
equation  expressing  stress  in  terms  of  strain  is  developed.  Index 
notation  closely  following  the  form  by  Fung  [39]  is  used  in  this 
derivation.  This  analysis  was  developed  in  terms  of  incremental 
stress  and  strain  for  a  generalized  plane  strain  formulation  by 
Adams  [11],  based  upon  work  by  Swedlow  [40],  The  terra  \  in  Eq. 

(A-l)  is  thus  developed  as  implied  in  Ei,  (A-2): 


e.  (p)  *  (  sklaki  }  s. . 

cij^  1  etoZMj.  >  siJ 


(A-2) 


where  the  dots  indicate  increments  of  stress  and  strain,  as  distin¬ 
guished  from  the  total  stress  or  strain  of  the  material,  t0  » 

(-j  s^js^j)  is  the  octahedral  shear  stress,  2Mf  is  the  tangent  modulus 
of  the  octahedral  shear  stress  -  octahedral  plastic  shear  strain  curve, 
and  represents  the  incremental  stress  tensor.  Adding  elastic 

incremental  strains  to  the  plastic  component  above  yields  the  overall 
expression  for  incremental  strain  [  1 1 ] , 


103 


t  .  . 

ij 


1  +  v  .  v  • 

-  —  a 


„  ,  siisk.l°kl 

O  .  .  +  - r - - 


ij  E  kk  ij  6t  2Mt 


(A- 3) 


!  obtain  an  equation  for  stress  in  terms  of  strain,  Eq .  (A-3)  must  be  in¬ 
verted.  Inverting  this  three-dimensional  form  follows  closely  the  procedure 
for  the  case  of  plane  strain  £  11  ]  .  The  first  step  is  to  express  the  second  and 
thira  terms  on  the  right-hand  side  of  Eq.  (A-3)  in  terms  of  strains  by 
multiplying  Eq.  (A-3)  by  the  Kronecker  delta  (6jj),  and  simplifying, 


ii  E  kk 


3  -  .  ,  Sj^SklCki 

—  a,  .  +  - 


(A-4) 


it  will  be  noted  that  slt  is  the  first  invariant  of  tiie  deviatoric 
itre^j  tens,  r,  wnich  is  identically  zero.  Thus,  Eq .  (A-4)  becomes 


ti 


(A-5) 


w:ucn  is  easily  inverted  t.>  obtain  an  expression  for  the  second  term: 


-)  •; 


it 


(A-6) 


1.  r  tiie  third  term,  F,q .  (A-3)  is  multiplied  by  the  deviatoric  stress 

•  .-ns.  r  and  simplified. 
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which  when  substituted  into  Eq.  (A-3)  along  with  Eq.  (A-6)  results 
in  the  desired  constitutive  equation: 
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Expanding  Eq. 
represented  by 
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A- 10)  in  matrix  form  results  in  the  more  useful  form 
Eq.  (6). 
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APPENDIX  B 


TWO-DIMENSIONAL  GENERALIZED  PLANE  STRAIN  FORMULATION 


The  displacement  components  must  be  in  the  following  form  for 
generalized  plane  strain,  as  previously  described: 

u  =  u(x,y) 

v  =  v(x,y)  (B-l) 

w  =  w(x,y)  +  C^z 

The  changes  involved  in  converting  from  plane  strain  to  generalized 
plane  strain  will  affect  the  strain-displacement  relation  and  sub¬ 
sequently  the  entire  stiffness  formulation.  Following  the  derivation 
by  Zienkiewicz  [21],  and  including  the  new  generalized  features  of 
plane  strain,  the  strain-displacement  behavior  of  a  triangular  element 
is  derived  below. 

The  form  for  the  desired  elemental  strain-displacement  relation¬ 
ship  wi] L  be: 

{F}  =  [N] (6)  (B-2) 

where  i F }  represents  displacements  at  any  point  within  the  element, 

(5}  represents  nodal  displacements  of  the  element,  and  the  [N]  are 
shape  functions,  the  general  functions  of  position.  Considering  Eqs. 
(B-l),  the  expression  in  Eq.  (B-2)  can  be  specified  in  terms  of 
generalized  plane  strain: 
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u(x,y) 
{F}  =(  v(x,y) 

w(x,y,  z) 


( B—  3 ) 


..  -  £ 
1  R 


( B— 4) 


where  the  displacement  subscript  represents  a  certain  node. 
To  solve  for  the  shape  functions,  a  linear  polynomial  is  chosen  for 
each  displacement  equation.  Thus,  the  "constant  strain"  element  is 


created,  as  defined  by 


u  =  Ji+a2X  +  ajy 


v  =  +  oi  cx  +  agy 


(  B—  5 ) 


w  =  ay  +  agX  +  agy  +  ^lO2 


Applying  nodal  conditions  to  the  expressions  for  u  and  v  follows  the 
work  of  Zienkiewicz  [21]  exactly,  resulting  in  representations  of 
displacements  in  terms  of  shape  functions  and  nodal  displacements: 


S'  =  (a.  +  b  x  +  c.y)/2A 


(B-6) 


where 


a.  =  X .  Y  -  X.Y. 
l  j  k  kj 


b.  =  Y.  -  Y. 
i  J  k 


c.  =  X.  -  X. 
1  x  j 


Shape  functions  N  and  follow  the  same  pattern.  It  will  be  noted 
that  the  shape  functions  are  functions  of  nodal  coordinates  and  ele¬ 


ment  area.  Expressions  for  u  and  v  become: 
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u  =  N.ui  +  NjUj  +  Nk\ 

v  =  N.v.  +N.v.  +N,  v, 
i  i  11  k  k 


(B-7) 


The  nodal  conditions  for  w  displacements  require  that  a  normal  load 
in  the  z-direction  induce  a  constant  displacement,  w  ,  regardless  of 
x  and  y  coordinate  positions.  Assuming  the  model  to  have  a  thickness 
of  unity, 

w  =  w.  @  z  =  0 

1  (B-8) 

w=w  +  w.  @  Z  =  1 
n  1 

Applying  these  conditions  to  the  third  of  Eqs .  ( B—  5) : 

w.  =07  +  ctgx  +  o9y 

1  ( B— 9 ) 

+  w  =  07  +  C19X  +  09V  +Oiq(1) 

Subtracting  the  first  of  Eqs.  (B-9)  from  the  second  results  in 

w  =  oin  (3-10) 

n 

Substituting  back  into  the  third  of  Eqs.  (B-5) : 

w  =  07  +  a gx  +  a9y  +  w^z  ( B— 11) 

Rearranging, 

w  -  w^z  =  07  +  a8x  +  a9y  (B-12) 

Solving  (B-12)  for  the  shape  functions  is  exactly  like  that  used  to  ob¬ 
tain  the  solutions  for  u  and  v,  resulting  in  the  same  form  of  expressions 
for  the  shape  functions 

w  -  w  z  =  N.w.  +  N.w.  +  N,  w,  (B-13) 

n  11  j  j  k  k 

Rearranging  gives  the  final  expression  for  w  in  terms  of  the  element 

shapes  and  nodal  displacements, 

w  =  N.w,  +  N.W.+N,  w.  +  wz  (5-14) 

1  i  11  k  k  n 


) 
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To  express  the  element  strains  in  terms  of  nodal  displacements,  the 
strain-displacement  relations  of  a  continuum  must  be  followed, 


3u 

y  ~ 

3u 

+ 

3v 

dX 

xy 

9y 

3x 

3v 

9u 

+ 

3w 

3v 

rxz 

5  z 

3x 

3w 

JV 

+ 

3w 

aZ 

v  z 

jZ 

3y 

Using  Eqs.  (B-7)  and  ( B— 14)  in  Eqs.  (B-15)  results  in  the  following 


Combined  into  matrix  form,  the  strain-displacement  relations  for 
generalized  plane  strain  are  defined  by 
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APPENDIX  C 


COMPUTER  PROGRAM 
FLOW  CHART 


MAIN 


An  administrative  program  which  calls  the  working 
subroutines  in  the  proper  order.  Reads  administra¬ 
tive  data,  i.e.,  titles,  fiber  volume,  scale  factors 
for  plocting,  number  of  loads,  output  flags,  initial 
temperature,  and  moisture  content.  Loop  on  load  in¬ 
crements  . 

GDATA 

Reads  mesh  data  for  the  finite  element  model, 
establishes  node  point  coordinates  and  connectivity, 
boundary  conditions  and  expected  bandwidth  of  the 
global  stiffness  matrix.  Builds  the  boundary 
condition  vector  and  indicating  vector  for  special 
elimination  operations  in  forming  the  stiffness  matrix. 

MESH 

Plots  the  finite  element  grid. 


FORMK 

Administrates  the  formation  of  the  global  stiffness 
matrix  by  calling  respective  subroutines  that  form 
elemental  stiffness  matrices.  Forms  global  stiff¬ 
ness  matrix  in  upper  triangular  form  and  simul¬ 
taneously  eliminates  rows  and  columns  where 


boundarv  conditions  are  involved. 
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Calculates  material  pro- 

1  Calculates  plasticitv 

perties  as  a  tine  t  ion  of 

j  constants. 

moisture  and  temperature. 

Octahedral  values  are 
available  lor  inelastic 
behavior  and  va. lure 
ana  lvsis . 

■  SOLVE 

Solves  simultaneous  equations  using  specialized 
Gaussian  .•! Im-nat ion  for  a  triangularined  matrix 
i  with  special  Sranca  [  1 0 ]  boundarv  condition 

I 

columns.  The  resulting  displacements  in  condensed 
form  are  expanded  to  the  original  size  hir  all 
nodes  for  subroutine  STRESS. 


STRESS 


Administrates  calculations  of  stress  and  failure 
criteria.  Calculates  and  outputs  incremental 
and  total  stresses  and  strains.  Checks  for 
unloading  elements.  Administrates  calculations 
of  principal  stresses  and  interface  stresses. 
Checks  for  hydrostatic  failure.  Return  to  MAIN 
for  next  load  increment. 


CMSTR 


Calculates  elemental  stresses  and  adds  hvgro- 
thermal  stresses. 


OCTA 


Checks  for  octahedral  failure  or  fiber  failure 
and  sets  failure  flags  for  STRESS.  Calculates 
failure  ratio  if  applied  stresses  are  higher 
than  failure  stresses. 


OCSTN 


Calculates  element  octahedral  shear  strains. 


PRINCE 

Calculates 

tensor . 

principal  values  for  stress  or  strain 
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